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ABSTRACT 


Multigrid  methods  were  developed  to  solve  partial  differential  equations. 
Research  has  shown  that  these  methods  are  applicable  to  a  broader  range  of 
problems.  This  thesis  investigates  the  application  of  multigrid  techniques  to 
minimal  cost  flow  problems,  specifically  the  long  transportation  problem. 

This  research  shows  that  multigrid  techniques  can  be  successfully  applied 
to  large-scale  long  transportation  problems  posed  on  a  three-dimensional,  regular 
grid  in  cost  space.  A  V-cycle  algorithm  is  developed  for  the  long  transportation 
problem.  Analogies  to  the  multigrid  components  of  restriction,  interpolation  and 
relaxation  are  detailed.  Performance  of  the  algorithm  is  discussed,  and 
computational  cost  is  analyzed. 

Future  research  is  likely  to  include  the  development  of  more  sophisticated 
restriction  and  interpolation  schemes  to  provide  integer-valued  flows,  and  the 
development  of  a  method  to  map  an  irregularly  spaced  problem  to  a  regular  grid, 
and  to  map  the  regular  grid  solution  back  to  the  original  problem  domain. 
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THESIS  DISCLAIMER 


The  reader  is  cautioned  that  the  computer  program  developed  during  this 
research  has  not  been  tested  on  all  cases  of  interest.  While  every  effort  has  been 
made  to  ensure  that  there  are  no  logic  or  computational  errors,  the  program 
cannot  be  considered  verified  as  yet.  Any  application  of  these  programs  without 
additional  verification  is  at  the  risk  of  the  user. 
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I.  INTRODUCTION 


A.  THE  TRANSPORTATION  PROBLEM 

The  generic  transportation  problem  contains  m  origins  and  n  destinations. 
Origins  and  destinations  are  also  referred  to  as  nodes  throughout  this  thesis.  Each 
origin  i  contains  a  supply  of  s,  units  of  a  particular  commodity.  Each  destination 
i  has  a  demand  for  d,  units  of  the  same  commodity.  It  is  assumed  that  total 
supply  equals  total  demand  and  that  s,,d^-  The  problem  is  to  find  a  feasible 
shipping  pattern  that  minimizes  total  cost.  The  problem  is  represented 
mathematically  as  follows; 


Minimize  2=c„a:„  +  ...  +  +  ... 


Subject  to:  +  ...  +  x,^ 


X  ,  ■*•...+  X  =  s 

ml  rrm  m 


+  x,,  +x 


m! 


X 


In 


♦  X. 


2n 


+  X 


X..  >  0  i=l,2,...,m, 


1 


In  general,  a  transportarion  problem  with  m  supply  nodes  and  n  demand  nodes 
contains  mn  variables  and  m+n  constraint  equations. 

A  transportation  problem  where  m«n  is  referred  to  as  the  long  transportation 
problem. 

B.  MULTIGRID  METHODS 

Multigrid  methods  evolved  out  of  an  effort  to  improve  iterative  solution 
methods  for  boundary  value  problems.  The  continuous  boundary  value  problem 
is  transformed  into  a  discrete  problem  by  choosing  a  set  of  grid  points  in  the 
domain  for  the  problem.  At  each  grid  point,  an  approximation  to  the  differential 
operator  is  formed  using  the  difference  between  values  of  the  unknown  function 
at  neighboring  grid  points.  This  leads  to  a  system  of  linear  equations  relating  the 
values  of  the  unknown  function  at  the  grid  points.  A  one-dimensional  problem 
is  used  in  this  thesis  to  illustrate  multigrid  methods  for  ease  of  demonstration, 
although  multigrid  techniques  are  more  commonly  applied  to  higher  dimensional 
problems. 

C  PREVIOUS  RESEARCH 

Aggregation/disaggregation  is  a  method  applied  to  large  scale  transportation 
problems.  It  involves  consolidating  several  neighboring  origins  or  destinations  to 
form  a  smaller  problem  that  is  easy  to  solve,  and  provides  an  initial  guess  for  the 
original  problem.  This  method  aggregates  demand  nodes  until  a  small  problem, 
which  is  easy  to  solve,  is  obtained.  The  solution  is  then  disaggregated  to  provide 
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an  initial  solution  or  a  bound  for  the  original  problem.  Kaminsky  (1989) 
developed  a  multilevel  algorithm  that  uses  an  approach  similar  to  the 
aggregation/disaggregation  approach.  Kaminsky  extended  this  approach  from 
two  levels  to  multiple  levels.  He  required  that  the  supply  and  demand  nodes 
occupy  a  physical  location  in  space  and  that  shipping  costs  be  directly  related  to 
the  physical  distance  between  a  supply  node  and  a  demand  node.  He  then 
aggregated  nodes  together  that  were  physically  close. 

Cavanaugh  (1992)  relaxed  tnis  geometric  restriction  so  that  shipping  costs 
may  have  no  correlation  to  the  physical  distance  between  the  supply  and  demand 
nodes.  He  then  mapped  the  demand  nodes  to  what  he  called  cost  space,  where  the 
shipping  cost  from  each  supply  node  determines  its  position.  For  example,  cost 
space  for  the  mxn  problem  is  the  m-dimensional  space  in  which  each  coordinate 
axis  represents  the  shipping  costs  from  each  of  the  m  supply  nodes.  Each  of  the 
n  demand  nodes  are  then  placed  in  cost  space  at  the  point  whose  coordinates  are 
the  shipping  costs  from  the  supply  nodes  to  it  (Cavanaugh,  1992).  This  generally 
results  in  a  transportation  problem  posed  on  an  irregular  grid  in  cost  space.  The 
idea  is  to  aggregate  the  demand  nodes  together  that  have  similar  costs. 

This  thesis  uses  the  idea  of  cost  space  and  assumes  that  the  problem  is 
posed  on  a  regular  grid. 
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D.  THESIS  OVERVIEW 


The  algorithm  developed  in  this  thesis  is  not  intended  to  replace  or  compete 
with  any  of  the  existing  solution  methods.  The  transportation  problem  is  being 
used  as  a  test  bed,  and  an  effective  method  would,  after  development,  be  applied 
to  a  class  of  more  intractable  problems. 

Chapter  II  provides  backgroimd  information  on  the  transportation  problem 
and  outlines  traditional  solution  methods.  Chapter  III  provides  an  introduction 
to  multigrid  techniques  and  describes  the  key  elements  in  multigrid  methods. 
Chapter  IV  describes  the  analogies  to  the  key  elements  in  multigrid  methods 
developed  in  this  thesis.  Chapter  V  describes  the  entire  algorithm  and  contains 
the  results  of  the  numerical  experimentation.  It  also  describes  the  storage 
requirements  and  provides  a  complexity  argument  of  the  algorithm.  Chapter  VI 
contains  conclusions  and  recommendations  for  future  research. 
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II.  INTRODUCTION  TO  LINEAR  PROGRAMMING 


A.  GENERAL  LINEAR  PROGRAMMING 

Linear  programming  is  concerned  with  the  optimization  of  a  linear  function 
with  n  variables  subject  to  m  linear  constraints.  The  general  form  of  the  linear 
programming  problem  is  as  follows: 

Minimize  Z"  cjr  (D 

;  ) 

Subject  to;  i=l,2,...,m  ^2) 

X.  >  0  (3) 

The  objective  function  is  represented  by  the  summation  in  (1)  and  is  usually 
denoted  by  2.  The  term  optimization  refers  to  the  minimization  or  maximization 
of  the  objective  function.  A  minimization  problem  can  be  easily  transformed  into 
a  maximization  problem  and  vice  versa  by  simply  taking  the  negative  of  the 
objective  function  as  follows:  MaxXjCjXj=  -MinLj(-CjX|).  The  cost  coefficients  are 
CiX^2f...,c„  and  the  decision  variables  are  The  constraint  inequalities  are 

represented  by  (2).  The  coefficients  a^  for  i=l,...,m  and  ;=I,...,n  are  referred  to  as 
the  technological  coefficients.  Inequality  (3)  is  the  nonnegativity  constraint.  The 
inequality  constraints  can  be  easily  transformed  into  equalities  by  adding  a  slack 
or  surplus  variable  0.  Hence  Za^,  S  bf  can  be  rewritten  as  la^i  -  =  bf, 


likewise  <  b,  can  be  rewritten  as  +  x,*,  =  b,.  A  set  of  values  for 
satisfying  all  the  constraints  is  referred  to  as  a  feasible  point,  feasible  vector,  or  a 
basic  solution,  although  it  may  not  yield  an  optimal  value  for  the  objective 
function.  The  set  of  all  such  points  is  known  as  the  feasible  region  (see  Figure  1). 
If  the  solution  contains  fewer  than  m+n-1  positive  solution  variables  then  the 
solution  is  said  to  be  degenerate. 


Figure  1.  Illustration  of  Feasible  Region,  the  Region  in  Which  All  Constraint 
Inequalities  are  Simultaneously  Satished 
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The  problem  can  be  presented  in  matrix  notation  as  follows: 


Minimize  cx 
Subject  to:  Ax  =b 

X  ^  0 


Note  that  the  inequality  constraints  have  been  transformed  into  equality 
constraints;  this  is  standard  practice,  since  systems  of  equations  are  easier  to  solve 
than  are  systems  of  inequalities. 

B.  MINIMAL  COST  FLOW  PROBLEMS 

A  graph  consists  of  a  set  V  of  m  nodes  (or  vertices)  and  a  set  £  of  edges,  which 
are  unordered  pairs  of  elements  from  V.  A  directed  graph,  or  digraph,  D,  consists 
of  a  set  V  of  m  nodes  and  a  set  A  of  arcs  which  are  ordered  pairs  of  elements 
from  V.  The  nodes  are  re,'  v  esented  by  points  and  there  is  a  directed  line  from  i 
to  ;  if  and  only  if  i,je  V  and  U,j)e  A.  A  network  is  a  directed  graph  in  which  every 
node  i  has  a  number  associated  with  it,  where  i=l,...,m,  and  each  arc  has  a 
capacity  for  flow  that  is  greater  than  or  equal  to  zero  (see  Figure  2).  The  number 


7 


r 


b,  represents  either  the  supply  of  some  commodity  present  at  the  node  or  the 
demand  of  the  commodity,  that  is,  the  amount  required  at  the  node. 


Figure  2.  Graphical  Representation  of  a  Mimimal  Cost  Flow  Network 


The  minimal  cost  flaw  problem  is  to  ship  the  available  supply  through  the  network 
to  satisfy  the  demand  at  the  cheapest  possible  cost.  Let  m  represent  the  total 
number  of  nodes  in  the  network.  Let  represent  the  flow  along  the  arc  from 
node  i  to  node  /.  Let  represent  the  cost  to  ship  one  unit  of  flow  along  the  arc 
(i,j).  Both  and  c^j  are  assumed  to  be  nonnegative.  If  b,<0,  then  b,  represents  the 
demand  at  node  i,  while  if  b,>0,  then  b,  represents  the  supply  at  node  i.  The 
underlying  assumption  in  this  problem  is  that  total  supply  equals  total  demand, 
therefore  'LJb-O.  If  this  is  not  the  case  then  a  dummy  node  can  be  added  to  the 
network  so  that  b^j=-'Z,bf,  and  then  arcs  with  zero  cost  from  each  supply  node  to 
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the  dummy  node  b^i  would  also  be  added.  The  general  minimal  cost  flow 
problem  can  be  written  mathematically  as  follows: 

Minimize  ZXr  jc 

I  I  1  ^ 

Subject  to:  ^  Vi  ^4) 

X.,  >  0, 

1 

Equation  (4)  is  referred  to  as  the  conservation  of  flow  equation.  represents  the 
total  flow  out  of  node  i,  and  l^Xj,  represents  the  total  flow  into  node  i.  Because  the 
network  is  actually  a  directed  graph,  solution  techniques  can  take  advantage  of 
this  special  structure  to  solve  these  problems  much  faster  than  general  linear 
programming  problems. 

C  THE  TRANSPORTATION  PROBLEM 

The  transportation  problem  is  the  simplest  case  of  the  minimal  cost  flow 
problem.  It  can  be  represented  by  a  complete  bipartite  digraph  (see  Figure  3).  This 
digraph  is  bipartite  because  the  set  of  nodes  can  be  partitioned  into  two  disjoint 
sets,  0={origin  or  supply  nodes}  and  D={destination  or  demand  nodes},  and  all 
arcs  are  directed  from  O  to  D.  The  digraph  is  complete  because  there  is  an  arc 
connecting  every  node  in  O  to  every  node  in  D.  If  an  actual  problem  does  not 
contain  an  arc  connecting  some  node  ie  O  to  some  node  ;e  D,  then  an  artificial  arc 
is  added  to  the  graph  and  assigned  high  cost  so  that  x^  becomes  an  artificial 
variable. 
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origins 


Destinations 


Figure  3.  Graphical  Representation  of  the  Transportation  Problem 

The  assumption  that  total  supply  equals  total  demand  guarantees  the 
existence  of  a  feasible  solution  to  the  transportation  problem,  provided  that  no 
constraints  are  placed  on  the  carrying  capacity  of  the  arcs.  For  example,  it  can  be 
shown  that  x~(s4j)ld  is  a  feasible  solution,  for  all  where 

Also,  a  bound  exists  for  each  decision  variable,  such  that, 
0^^^minimum{s^,d^}.  It  is  well-known  that  a  bounded  linear  program  with  a 
feasible  solution  has  an  optimal  solution  (Bazaraa,  1990).  Some  traditional  solution 
methods  are  presented  next. 


D.  THE  SIMPLEX  METHOD 

Whenever  feasible  solutions  to  a  linear  programming  problem  exist,  the 
region  formed  by  the  constraint  equations  is  called  the  feasible  region  (see  Figure 
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1).  It  can  be  shown  that  this  region  is  also  convex,  which  means  that  each  pair  of 
points  in  the  region  can  be  joined  by  a  line  segment  that  lies  within  the  region. 
The  boundaries  of  this  region  are  lines  or  planes  and  there  are  comers  where 
these  lines  or  planes  intersect.  Points  on  these  corners  are  referred  to  as  extreme 
points.  Two  extreme  points  are  adjacent  if  the  line  segment  joining  them  is  an  edge 
of  the  feasible  region.  Any  point  in  the  region  can  be  written  as  a  convex 
combination  of  the  extreme  points,  which  leads  to  the  property  that  an  optimal 
solution  can  be  found  at  an  extreme  point. 

Informally,  the  simplex  method  can  be  summarized  as  follows: 

1.  Begin  at  an  initial  extreme  point. 

2.  Perform  an  optimality  test  (Is  this  extreme  point  at  least  as  good  as  the 
adjacent  extreme  point(s)?).  If  yes,  stop;  optimality  has  been  achieved. 

3.  If  not,  move  to  an  adjacent  extreme  point  whose  objective  function  value 
is  at  least  as  good  as  the  present  extreme  point.  Repeat  step  2. 

The  simplex  method,  as  presented  above,  is  an  algorithm  used  to  solve  exactly  a 

linear  programming  problem  in  a  finite  number  of  steps  or  reveal  that  the 

solution  is  unbounded,  that  is,  that  the  objective  function  has  a  value  of  -«»  and 

hence  no  optimal  solution  exists.  This  method  breaks  a  difficult  problem  into  a 

series  of  easy  problems.  Some  notation  must  be  introduced  before  the  simplex 

method  can  be  formally  presented. 

Consider  the  system  Ax=b,  x^.  Let  A  be  an  mxn  matrix,  where  m^,  and 
b  a  vector  of  length  n,  and  suppose  that  the  rank  of  the  augmented  matrix  [A,b] 
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sarisfies  rank  ([A,b])=rank(A)=m.  The  matrix  A  can  be  partitioned  into  A=(B  NJ, 
where  B  is  an  mxm  invertible  matrix  and  N  is  an  mx(n-m)  matrix.  A  solution  to 
the  system  is  x^=Ixb  x^l,  where  Xb=  B  'b  and  X|^=0.  This  is  referred  to  as  a  bask 
feasible  solution.  B  is  the  bask  matrix  and  N  is  the  nonbask  matrix.  The  columns  of 
B  are  also  referred  to  as  the  basis.  The  vector  x,  contains  values  of  the  bask 
variables  and  x^  contains  values  of  the  nonbask  variables.  If  Xg>0  then  the  solution 
is  nondegenerate,  while  if  at  least  one  component  of  x,  equals  zero  then  the 
solution  is  degenerate.  The  cost  vector  is  partitioned  so  that  c=[cg  0^1,  where  c, 
contains  the  basic  variable  costs  and  contains  the  nonbasic  variable  costs.  The 
dual  variables  are  contained  in  the  row  vector  u,  where  u  is  given  by  u=c,B  '.  The 
dual  variables  are  used  to  compute  the  reduced  costs  denoted  by  (2y-cp=uaj-Cy, 
where  is  the  jth  column  of  the  matrix  A  and  is  given  by  z-c^B  ^  for  each 
nonbasic  variable.  The  reduced  costs  are  actually  the  cost  coefficients  of  the 
nonbasic  variables.  The  simplex  algorithm  for  a  minimization  problem  is  as 
follows: 

1.  Start  with  a  basic  feasible  solution  with  basis  B.  (There  are  many  different 
procedures  available  for  finding  an  initial  basis.) 

2.  Solve  the  system  Bxg^b,  which  has  the  unique  solution  x,=B  ’b.  Let  Xn=0 
and  z=CbXb. 

3.  Solve  the  system  uB=Cb/  which  has  the  unique  solution  u=CbB  '.  Calculate 
Zy-c-uaj-c^,  for  all  nonbasic  variables.  Let  z^-c,=  maximum(Zy-Cy)  for  all  je R 
where  R  is  the  current  set  of  indices  associated  with  nonbasic  variables.  If 
Zj-Ci^O  then  stop,  the  solution  is  optimal.  Otherwise  x^  is  the  entering 
variable. 
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4.  Solve  the  system  By,,=ak,  with  unique  solution  y,^=B'*ak.  If  that  is,  if 
all  components  of  the  vector  y,,  are  negative,  then  the  optimal  solution  is 
unbounded.  Otherwise  continue. 

5.  The  variable  x^,  enters  the  basis,  while  leaves  the  basis.  The  index  r  of 
the  variable  Xg,  which  leaves  the  basis  is  determined  by  the  minimum  ratio 
test.  If  the  columns  of  B  are  the  basis  then  a^  enters  the  basis  while  x,, 
becomes  one  of  the  basic  variables.  The  minimum  ratio  test  is  as  follows; 


Br 


ykr 


=  minimum 
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6.  Update  the  basis  B,  where  a^  replaces  a^,,  update  the  set  R,  and  repeat  step 

2. 

A  transportation  problem  with  25  supply  nodes  and  1000  demand  nodes 
consists  of  25000  variables  and  1025  constraints.  This  type  of  problem  cannot  be 
solved  efficiently  or  quickly  using  the  simplex  method,  even  with  the  aid  of  a 
large  scale  computer.  The  network  simplex  method  is  suited  to  the  task  because  it 
exploits  the  unique  struchire  of  the  graphical  representation  of  the  transportation 
problem. 

E.  THE  NETWORK  SIMPLEX  ALGORITHM 

Several  terms  must  be  defined  before  the  network  simplex  method  can  be 
discussed.  A  digraph  D  contains  a  path  from  s  to  f  if  an  ordered  set  of  the  arcs  in 
A  begin  at  s  and  end  at  t,  where  s,f€V.  That  is,  the  arcs  (Wo,Wj),  (Wj,W2), 
(w2,W3),...,(w„.j,w„)  where  Wo=s  and  w„=t,  and  w^  are  all  in  A,  for;=0,2,2,...«.  A  path 
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that  begins  and  ends  at  the  same  node  is  a  cycle.  A  digraph  is  connected,  if 
between  every  pair  of  vertices  s  and  t  there  is  a  path.  A  tree  is  a  connected 
digraph  with  no  cycles.  A  spanning  tree  in  D  is  a  tree  containing  every  node  in  D. 
A  rooted  tree  has  one  node  identified  as  the  root  such  that  there  exists  a  unique 
path  from  the  root  to  every  node  in  D. 

If  the  solution  to  a  minimal  cost  flow  problem  is  examined  graphically  it 
corresponds  to  a  spanning  tree  in  the  network  (Bazaraa,  1990).  The  network 
simplex  method  exploits  this  fact  and  is  a  more  efficient  solution  method  than  the 
simplex  method.  The  network  simplex  algorithm  is  as  follows: 

1.  Find  an  initial  basic  feasible  solution  represented  by  a  rooted  spanning 
tree. 

2.  Compute  the  basic  flows,  x^,  and  the  dual  variables 

3.  For  each  nonbasic  arc  compute  the  reduced  cost 

4.  If  then  stop,  as  the  solution  is  optimal.  Otherwise  select  a 

nonbasic  arc  (i,;)  with  a  positive  reduced  cost  to  enter  the  basis  (append  to 
the  current  basic  tree). 

5.  Divert  flow  to  the  entering  arc  until  flow  along  another  arc  is  reduced  to 
zero.  Remove  the  arc  with  zero  flow  from  the  tree,  thereby  aeating  a  new 
basic  spanning  tree.  Go  to  step  2. 

This  method  can  be  applied  directly  to  the  actual  network  representation  of  a 
small  problem  and  it  can  be  programmed  to  solve  larger  and  more  complicated 
problems.  Recall  that  the  transportation  problem  is  the  simplest  case  of  a  minimal 
cost  flow  problem.  Large  scale  transportation  problems  are  normally  solved  using 
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the  network  simplex  method,  but  occasionally  a  method  known  as 
aggregation/disaggregarion  is  utilized. 

F.  AGGREGATION/DISAGGREGATION 

Aggregation/  disaggregation  (referred  to  as  aggregation)  is  a  method  used 
to  solve  large  scale  transportation  problems.  It  involves  partitioning  the  variables 
in  a  problem  and  replacing  each  partition  with  a  single  variable.  The  same 
technique  is  applied  to  the  constraints.  The  resulting  problem  is  referred  to  as  an 
aggregate  problem.  The  method  involves  consolidating  several  "neighboring" 
origins  or  destinations  to  form  a  much  smaller  problem  that  is  easier  to  solve  and 
that  will  provide  insight  into  the  solution  to  the  larger  problem.  The  term 
"neighboring"  is  intentionally  loosely  defined  so  that  there  are  not  explicit 
restrictions  on  the  aggregation  method  (Balas,  1963). 

Aggregation  provides  an  approximate  solution  to  the  original  problem  while 
considering  only  a  small  portion  of  the  problem  data.  To  solve  the  original 
problem  directly  using  network  simplex  with  m  origins  and  n  destinations 
requires  mn  evaluations  of  the  reduced  costs,  while  if  the  aggregation  method  is 
used  the  number  of  computations  decreases  significantly,  depending  on  how 
many  aggregate  variables  are  used.  For  instance,  if  a  problem  containing  m=350 
origins  and  n=500  destinations  is  solved  directly  it  w:’l  contain  mn=175,000 
variables  and  850  constraint  equations.  If  the  aggregation  method  is  used  and 
aggregation  is  five  by  five  (five  neighboring  origins  become  one  origin  and  five 
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neighboring  destinarions  become  one  desrination),  then  there  are  M=70 
aggregated  origins  and  iV=100  aggregated  destinations,  resulting  in  MN=7000 
variables  and  170  constraint  equations.  The  gain  in  computational  savings 
increases  with  the  size  of  the  problem  and  the  concentration  of  the  origins  and 
destinations.  This  method  is  not  widely  used,  because  the  accuracy  of  a  solution 
is  often  traded  for  computational  savings,  and  an  optimal  solution  is  therefore  not 
guaranteed.  It  must  be  emphasized  that  while  the  nontraditional  solution 
approach  to  the  long  transportation  problem  in  this  thesis  is  related  to  the 
aggregation  method,  there  are  significant  differences  between  them. 
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III.  AN  INTRODUCTION  TO  MULTILEVEL  TECHNIQUES 

A.  THE  MODEL  PROBLEM 

Consider  the  boundary  value  problem  which  describes  the  steady  state 
distribution  of  temperature  in  a  long  uniform  rod.  This  may  be  represented  by  the 
second  order  differential  equation  -u''(x)+au(x)=f(x),  where  0<x<l,  aW  (5)  and 
subject  to  the  boundary  conditions  u(0)=u(l)=0  (Briggs,  1987).  This  equation  will 
be  referred  to  as  the  model  problem. 

Occasionally,  the  model  problem  can  be  solved  analytically,  but  it  is  more 
common  to  solve  it  numerically.  Multigrid  methods  are  applied  to  numerical 
solution  techniques.  One  such  numerical  method  of  solution  is  the  finite  difference 
method.  The  domain  0:^<1  is  partitioned  into  N  subintervals  with  spacing  h  = 
1/N,  and  each  point  is  labeled  Xj,  where  Xy=  jh.  This  is  the  finest  grid  and  is 
denoted  by  D**  (see  Figure  4). 


Figure  4.  One-dimensional  grid  on  the  interval  0<x<l.  Grid  spacing  is  h=l/N. 
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At  each  grid  point  the  differential  equation  is  replaced  with  a  finite  difference 
equation.  Let  u(x^  represent  the  exact  solution  and  v(x^  represent  the 
approximation  to  the  exact  solution,  where  For  ease  of  notation  let 

represent  u(Xj)  and  represent  v(Xj),  and  let  v=(t;,,  Vj,  We  seek  to  find  v 

such  that  the  components  of  v  satisfy  N-J  linear  equations,  each  having  the  form 
(•Vi.^+2vi-Vi^i)/h^+aVj=f(Xj),  where  The  boundary  conditions  yield  Vo=Vtf=0. 

This  discretization  may  be  represented  in  matrix  notation  as 


’2-kt/i^  -1 

-1  2-kj/i^  -1 

f^ 

-1 

• 

- 

-1  2^h\ 

or  simply  as  Av  =  f  (Briggs,  1987).  Note  that  the  matrix  A  is  tridiagonal,  sparse 
and  symmetric  positive  definite  with  dimensions  (N-l)x(N-l).  This  system  could 
be  solved  directly  using  Gaussian  elimination,  but  this  method  is  computationally 
expensive.  Indirect  methods  such  as  iteration  can  also  be  used  to  solve  systems  of 
this  type.  A  special  case  of  iteration  is  relaxation.  Jacobi  and  Gauss-Seidel  are  two 
examples  of  relaxation  methods.  Relaxation  begins  with  an  initial  guess  at  the 
solution,  and  then  tries  to  improve  the  current  approximation  through  a  series  of 
iterations.  Ideally  this  sequence  of  approximations  will  converge  to  the  exact 
solution.  The  rate  of  convergence  of  basic  iterative  methods  for  certain  problems 
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can  be  improved  significantly  using  multigrid  techniques.  A  discussion  of  basic 
iterative  methods  is  required  before  multigrid  techniques  can  be  presented. 


B.  BASIC  ITERATIVE  METHODS 

Let  Au  =  f  denote  the  linear  system  in  equation  (5).  The  exact  solution,  u, 
to  Au=f  is  unknown  but  the  approximation  v  is  known.  There  are  two  important 
measures  of  the  quality  of  the  approximation  v.  The  first  is  the  error,  denoted  by 

e=u-v.  The  size  of  the  error  can  be  measured  by  any  vector  norm  such  as 

I  lei  Since  the  error  cannot  be  found  if  the  exact  solution  is 

unknown,  consider  another  important  measure,  the  residual.  The  residual  r  is 

denoted  by  r^f-Av;  it  represents  the  amount  by  which  the  approximation  fails  to 

satisfy  the  equation.  The  residual  can  be  measured  by  a  vector  norm  in  the  same 

manner  as  the  error.  Rewrite  rs=f-Av  as  Av^f-r  and  subtract  it  from  the  equation 

Au=f.  This  results  in  the  equation  Ae=r,  referred  to  as  the  residual  equation. 

Multigrid  methods  employ  iterative  improvement,  that  is,  the  residual  equation 

is  solved  (approximately)  for  e  and  a  new  approximation  is  computed  to  improve 

V  in  the  next  iteration  using  the  definition  u=v+e.  That  is,  v‘"*’'=v'"*+  C,  where  C 

is  an  approximation  to  e. 

Consider  the  Jacobi  relaxation  method  when  applied  to  the  model  problem, 
rewritten  as  -M^,+2uy-Wy+,=/iy^,  where  Uo=Ufj=d)  and  o=0.  Let  be  the 

current  approximation  and  v'"*”  be  the  updated  approximation.  The  Jacobi 
method  consists  of  solving  the  f  equation  for  the  f  unknown  u^,  holding  all  of 
the  other  unknowns  fixed.  When  this  has  been  done  for  all  N-2  equations,  we 
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have  performed  one  sweep  of  the  Jacobi  iteration.  That  is,  for  j=l  compute 

1  ^j  <N-l .  Ideally,  these  iteration  sweeps  are  computed  until 

convergence  to  the  solution  is  attained. 

Typically,  after  a  few  sweeps  the  convergence  stalls,  meaning  that  the 
answer  will  not  improve  with  more  iterations.  This  is  because  the  error  is  made 
up  of  oscillatory  and  smooth  components,  and  many  iterative  methods  eliminate 
the  oscillatory  components  of  the  error  while  having  very  little  effect  on  the 
smooth  components.  This  is  known  as  the  "smoothing  property."  The  problem  is 
now  how  to  reduce  the  smooth  components  of  the  error  that  stiU  exist.  Multigiid 
methods  evolved  from  attempts  to  solve  this  problem.  The  following  sections 
address  this  problem  in  greater  depth. 

C  FOURIER  MODES 

A  short  discussion  of  Fourier  modes  is  necessary  to  illustrate  a  key  idea  in 
multigrid  methods.  Consider  the  model  problem  -u  "(x)=0.  This  is  discretized  and 
yields  a  simple  linear  homogenous  system  Au  =  0.  The  unique  solution  is  u=0 
and  the  error  in  the  approximate  solution  v  is  -v.  The  linear  system  Au=0 
becomes  where  Uf,=Uf^  and  Let  the  initial  guess  (and 

hence  the  initial  error)  be  one  the  Fourier  modes,  that  is,  a  vector  w^,  whose  f' 
component  is  w^.  -sin|~j  when  and  it  is  between  3  and  N-3.  The  integer 
k  is  the  wavenumber  and  it  represents  the  number  of  half-sine  waves  which 

constitute  on  the  domain  of  the  problem.  Small  values  of  k  correspond  to  long 
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smooth  waves  and  large  values  of  k  correspond  to  highly  oscillatory  waves. 
Figure  5  displays  several  Fourier  modes  on  a  grid  with  N=22.  The  mode 
consists  of  k  half  sine  waves  and  has  wavelength  l=2/k  (recall  that  the  length  of 
the  interval  is  1).  The  k  =  N/2  mode  has  a  wavelength  of  4h  and  the  k  =  N-1  mode 
has  a  wavelength  of  2h.  Waves  with  wavenumbers  greater  than  N  or  wavelengths 
less  than  2h  cannot  be  represented  on  the  grid.  Through  a  phenomenon  known 
as  "aliasing",  a  wave  with  length  less  than  2h  actually  appears  with  a  wavelength 
greater  than  2h.  Wavenumbers  where  1  <k  <  N/2  are  low  frequency  or  smooth 
modes,  and  wavenumbers  where  N/2  <k<  N-1  are  high  frequency  or  oscillatory 
modes. 

For  the  cases  Wj  and  w,2,  ten  relaxation  sweeps  are  performed  on  a  grid  with 
N=64  (see  Figure  6).  After  ten  relaxation  sweeps,  the  high-frequency  modulation 
on  the  long  wave  is  nearly  eliminated,  but  the  smooth  component  remains.  The 
key  point  of  this  discussion  is  that  there  is  a  fast  rate  of  convergence  as  long  as 
the  error  has  high  frequency  components  and  convergence  seems  to  stall  when 
only  low  frequency  components  remain. 

This  supports  the  earlier  observation  that  for  oscillatory  waves  the  error  is 
damped  very  easily,  but  the  smooth  components  remain  essentially  unchanged 
and  so  convergence  stalls.  The  multigrid  approach  is  to  treat  the  smooth 
component  of  the  error  on  a  coarser  grid  where  it  will  appear  more  oscillatory 
and  can  be  damped  so  that  fast  convergence  to  the  solution  is  possible. 
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Figure  5.  Graphs  of  the  Fourier  modes  on  a  grid  with  N=I2,  for  wavenumbers 
jt=I,4,  9. 

D.  COMPONENTS  OF  MULTIGRID 


Suppose  a  relaxation  scheme  is  applied  until  only  the  smooth  error 
components  remain.  The  smooth  wave  on  the  grid  O'*  looks  more  oscillatory  on 
the  grid  This  is  because  the  mode  k  is  the  same,  but  it  is  the  out  of  N/2, 
rather  than  out  of  N.  Also  the  grid  points  on  are  the  even  numbered  points 


of  D'’.  Consider  the  Jt"*  mode  of  D**  evaluated  at  the  even-numbered  grid  points. 
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Figuie  6.  Initial  guess  (a)  Wj,  (b)  w,2  and  (c)  combination  of  W2  and  w,2>  The 
figures  on  the  left  show  the  approximation  before  iteration.  The  figure  on  the 
right  is  after  ten  iterations. 

Then  if  its  components  may  be  written  as 

iw*  -sin(.?^)  -sin(-^!L)  1^.^ 

^  N  N/2  2 

One  can  say  that,  under  restriction,  the  k*  mode  on  Q*'  becomes  the  k***  mode  on 
Therefore,  transferring  from  Gl*'  to  the  it"'  mode  becomes  more  oscillatory. 
This  is  true  for  wavenumbers  k,  where  l^k^J2  (see  Figure  7). 
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Figure  7.  A  wave  with  wavenumber  Jt=4  on  £1'’(N=12)  projected  onto  ii^(N=6). 


For  wavenumbers  k,  where  k>N/2,  the  oscillatory  waves  on  Q*'  are  aliased  to 
become  smooth  on  This  implies  that  the  high  frequency  error  must  be 
damped  before  moving  to  the  coarse  grid.  Twc  important  features  to  note  are: 

1.  Because  the  error  is  smooth  after  relaxation,  it  can  be  represented 
accurately  on 

2.  The  smooth  error  on  appears  more  oscillatory  on  where  it  is  more 
easily  eliminated,  since  relaxation  is  cheaper. 

This  leads  to  the  conclusion  that  an  improvement  to  iterative  methods  would  be 

to  transfer  the  smooth  error  components  to  the  coarse  grid,  where  the  error 

modes  appear  more  oscillatory,  and  to  then  perform  relaxation.  Recall  that  there 

exists  an  equation  for  the  error,  namely  e=u-v,  satisfying  Ae=r=f-Av.  Therefore 

the  residual  equation  can  be  used  to  relax  directly  on  the  error.  Relaxing  on  Au=f 
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with  an  initial  guess  v  is  equivalent  to  relaxing  on  the  residual  equation  Ae=r 
with  initial  guess  u-v. 

Nested  iteration  and  coarse  grid  correction  are  two  techniques  that  utilize  the 
coarse  grids.  Coarse  grid  correction  uses  the  coarse  grid  to  improve  the  current 
approximation.  Nested  iteration  uses  the  coarse  grid  to  improve  the  initial  guess. 

An  important  question  now  arises:  How  does  one  move  from  the  fine  grid 
n*'  to  a  coarser  grid  The  next  two  sections  address  this  question  of  intergrid 
transfers. 

1.  Interpolation 

The  method  used  to  transfer  information  from  the  coarse  grid  to  the 
fine  grid  is  interpolation,  also  referred  to  as  prolongation.  There  are  several 
interpolation  methods  to  choose  from,  but  only  linear  interpolation  will  be 

discussed  here.  is  the  linear  interpolation  operator.  Consider  only  coarse  to 
fine  grid  spacings  with  a  ratio  of  2:1.  is  a  linear  operator  that  transforms  coarse 

grid  vectors  into  fine  grid  vectors,  that  is  For  the  model  problem  this 

means 

Figure  8(a)  illustrates  this  transformation. 
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2.  Restriction 


The  method  used  to  transfer  information  from  the  fine  grid  to  the 
coarse  grid  is  restriction.  There  are  also  several  restriction  methods  available,  but 

only  full-weighting  will  be  presented  here.  The  restriction  operator  /“  transforms 


fine  grid  vectors  into  coarse  grid  vectors,  that  is  /“  v’'=d^.  For  the  model 
problem  this  means 


This  method  takes  a  weighted  average  of  neighboring  fine  grid  points.  Figure  8(b) 
illustrates  this  transformation. 

3.  Coarse  Grid  Correction 

Coarse  grid  correction  incorporates  the  use  of  the  residual  equation  to 
relax  on  the  error.  This  method  requires  relaxation  on  the  fine  grid  until 
convergence  stalls,  then  relaxes  on  the  residual  equation  on  the  coarse  grid  to 
obtain  an  approximation  to  the  error,  then  returns  to  the  fine  grid  to  correct  the 
initial  approximation.  Coarse  grid  correction  is  a  two  grid  process,  and  may  be 
written  as: 


1.  Relax  on  A*'u''=f*'  to  get  t^'. 

2.  Restrict  r“=Jf  (/'■-/iV). 

3.  Solve  A^^=r^  to  get  e^. 
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4.  Interpolate  and  correct  (Briggs,  1987). 


Figure  8.  (a)  Restriction  of  a  fine  grid  vector  to  a  coarser  grid,  (b)  Interpolation  of 
a  coarse  grid  vector  to  a  fine  grid. 


This  leads  to  the  question:  How  do  we  solve  for  the  coarse  grid  scheme? 

The  solution  is  to  recursively  apply  coarse  grid  correction.  This  leads  to  the  V- 
cycle. 
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4.  The  V-Cyde 

The  V-cyde  algorithm  applies  the  coarse  grid  correction  scheme 
recursively  until  the  coarsest  grid  is  reached,  and  a  direct  solution  to  the  residual 
equation  is  possible.  It  then  works  its  way  back  to  the  finest  grid.  Note  that  as 
this  correction  scheme  is  applied  recursively  to  solve  A^e^=r^  the  notation  can 
be  simplified  by  letting  1^=^*  and  e^=v^.  The  coarse  grid  correction  scheme 
applied  recursively  is  presented  next.  Let  Q  denote  the  number  of  grids,  where 
Q>h  and  let  L=2^’  so  that  the  coarsest  grid  has  sparing  Lh. 

Relax  on  u,  times  with  initial  guess  v\ 

Compute  /“  r^. 

Relax  on  A^^=P  u,  times  with  initial  guess  v^=0. 

Compute  f^=  r^. 

Relax  on  A^u**’=f^  u,  times  with  initial  guess  v^=0. 

Compute  r***. 

Solve  A'^u'^=P. 

Correct  v^^v^+  v®**. 

Relax  on  A^^=f^  Uj  times  with  initial  guess  v^. 

Correct  v^<-v^  +  /“  v^. 

Relax  on  A^u^=f^  Uj  times  with  initial  guess  v^. 

Correct  v^. 

Relax  on  A'Hi‘’=f’'  Uj  with  initial  guess  v**  (Briggs,  1987). 
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The  V-cycle  scheme,  where  v*'<-MV''(v*',f‘')  may  be  defined  recursively  to  give  a 
compact  description  as  follows: 

1.  Relax  o,  times  on  with  a  given  initial  guess  v**. 

2.  If  equals  the  coarsest  grid,  then  go  to  4. 

Else  P«-/“(f^-AV). 

v^i-0 

3.  Correct  v*'<-v^+/^v^. 

4.  Relax  u,  times  on  A^''=f''  with  initial  guess  (Briggs,  1987). 

The  V-cycle  is  one  of  many  multigrid  cycling  schemes.  It  is  depicted  in  Figure 
9(a). 

Recall  that  the  smoothing  property  is  effective  at  eliminating  the 
oscillatory  components  of  the  error  while  having  little  or  no  effect  on  the  smooth 
components.  We  use  coarse  grid  correction  to  eliminate  the  smooth  error 
components.  It  is  the  combination  of  these  techniques,  relaxation  for  oscillatory 
error  and  coarse  grid  correction  for  smooth  (./ror,  that  makes  multigrid  such  an 
effective  method. 

It  stands  to  reason  that  if  the  initial  guess  is  good,  then  the  iterative 
scheme  has  less  work  to  do  in  solving  the  problem.  A  useful  technique  is  to  first 
solve  the  problem  on  a  coarse  grid  to  obtain  a  good  initial  guess  on  the  fine  grid. 
This  method  is  called  nested  iteration. 
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5.  Nested  Iteration 


Nested  iteration  uses  the  coarsest  grid  to  obtain  an  initial  guess  to  the 
next  finer  grid  and  continues  this  process  unril  the  original  grid  is  reached.  Given 
the  original  problem  Au=f  posed  on  the  coarsest  grid,  nested  iteration  proceeds 
as  follows:  Solve  Au=f  on  the  coarsest  grid.  Interpolate  the  solution  to  the  next 
finer  grid  to  obtain  an  initial  guess.  ...(repeat  the  process)...;  solve  Au=f  on  Q*''and 
interpolate  to  obtain  an  initial  guess  on  solve  Au=f  on  and  interpolate  to 
obtain  an  initial  guess  on  and  solve  Au=f  on  ii''  to  obtain  the  solution  (Briggs, 
1987).  The  algorithm  that  combines  nested  iteration  with  the  V-cycle  is  known  as 
the  Full  Multigrid  (FMV)  V-cycle. 

6.  The  Full  Multigrid  V-Cycle 

An  efficient  means  of  obtaining  a  good  initial  guess  is  to  solve  the 
problem  first  on  a  coarser  grid  and  interpolate  that  solution  for  use  as  a  fine-grid 
initial  guess.  This  idea  is  nothing  more  than  nested  iteration.  The  next  question 
to  arise  is:  How  can  one  solve  the  coarse-grid  problem?  One  obvious  choice  is  to 
use  a  V-cycle.  This  leads  to  the  Full  Multigrid  (FMV)  scheme.  The  FMV  scheme 
where  v‘'<-FMV''(v*’,f*')  and  f^P',...,V',v^,...  are  set  equal  to  zero,  is  as  follows: 

1.  If  fl’’  is  the  coarsest  grid,  then  go  to  step  3. 

Else  F'<-/“(f'’-AV) 

FMV^(v^,P’). 

2.  Correct  v**  <-  v*'+/^v^. 
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3.  v^<r-  MV'’(v^f’)  \)o  times  (Briggs,  1987). 

The  FMV  is  depicted  in  Figure  9(b).  Full  multigrid  is  the  culmination  of  the  ideas 
presented  in  the  preceding  sections  in  this  chapter.  Multigrid  combines  each  of 
the  methods  just  described  into  an  extremely  efficient  algorithm. 


Figure  9.  (a)  The  V-Cycle.  (b)  The  FMV  V-Cycle.  The  respective  algorithms  visit 
the  various  grids  in  the  sequence  indicated. 


E.  MULTIGRID  METHODS  APPLIED  TO  THE  LONG  TRANSPORTATION 
PROBLEM 

One  of  the  objectives  of  this  thesis  is  to  develop  a  multigrid  algorithm  to 
solve  the  long  transportation  problem.  Analogies  to  the  components  of  multigrid 
must  be  found.  Recall  that  the  key  components  in  multigrid  are  interpolation. 
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restriction,  and  relaxation  which  are  combined  together  to  form  the  V-cycle.  One 
major  stumbling  block  is  in  drawing  an  analogy  between  some  component  of  the 
transportation  problem  and  the  residual  equation,  which  is  key  to  multigrid 
solution  methods.  The  following  chapter  presents  the  components  of  the 
algorithm  developed  during  this  thesis  research. 
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IV.  MULTIGRID  COMPONENTS  IN  THE 
LONG  TRANSPORTATION  PROBLEM 

A.  BACKGROUND 

As  stated  previously,  a  necessary  assumption  for  the  algorithm  developed 
here  is  that  the  problem  is  posed  on  a  regular  grid.  For  instance,  consider  a 
problem  with  m  supply  nodes  and  n  demand  nodes  posed  on  a  regular  m- 
dimensional  grid.  The  dimension  of  the  grid  is  determined  by  the  number  of 
supply  nodes,  hence  m  supply  nodes  implies  an  m-dimensional  regular  grid.  The 
point  identified  by  the  three-tuple  (i,j,k)  denotes  a  particular  demand  node.  The 
index  i  denotes  the  cost  to  ship  to  demand  node  (i,j,k)  from  supply  node  1;  the 
index  j  denotes  the  cost  to  ship  from  supply  node  2  and  the  index  k  denotes  the 
cost  to  ship  from  supply  node  3.  This  can  be  written  as  c„^u=i,  C2f^y=7/  and  C3,i^u=k. 
Therefore,  m  supply  nodes  correspond  to  an  m-tuple  identifying  each  demand 
node.  A  three  supply  node  problem  with  shipping  costs  mapped  to  a  regular  grid 
is  displayed  in  Figure  10. 

The  method  developed  here  takes  a  problem  posed  on  D'*  and  transforms 
it  into  a  problem  on  D**' where  it  is  easier  to  solve,  it  then  interpolates  the  answer 
back  to  C2*'.  Three  questions  immediately  come  to  mind.  The  first  is:  How  is  the 
problem  transferred  from  D**  to  The  second  is:  How  is  the  problem  solved 
on  12“*?  And  finally:  How  is  the  problem  transferred  from  £2^  to  back  to  £2**? 
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Figure  10.  A  Three  Supply  Node  Problem  Posed  on  a  Regular  Grid 

These  questions  are  all  answered  in  the  next  four  sections.  But  first  an  explanation 
of  the  grid  and  demand  node  relationship  is  necessary. 

The  finest  grid  in  this  algorithm  contains  n=(2'-ir  demand  nodes.  The 
integer  t  represents  the  number  of  grid  levels,  m  represents  the  dimension  or  the 
number  of  supply  nodes.  In  what  follows,  we  assume  m-3;  therefore,  the  most 
general  problem  considered  will  have  size  3xn.  The  root  of  the  t**'  grid  is  denoted 
by  (2^-1),  where  r=l,...t.  Table  1  illustrates  the  relationship  between  the  grid 


number,  the  root  of  the  grid  and  the  number  of  demand  nodes  for  up  to  five 
grids. 

TABLE  1.  RELATIONSHIP  BETWEEN  GRIDS  AND  DEMAND  NODES 


Gridt 

Dimension,  m 

Root,  2*-l 

Demand  nodes 
N=(2‘-l)“ 

1 

3 

1 

1 

2 

3 

3 

27 

3 

3 

7 

343 

4 

3 

15 

3375 

5 

3 

31 

29,791 

As  shown  above  the  number  of  grid  levels  is  directly  related  to  the  number  of 
demand  nodes  in  the  original  problem.  If  four  grids  are  chosen  then  r=l,...,4.  The 
finest  grid  corresponds  to  r=4  and  contains  a  3x3375  transportation  problem.  The 
coarsest  grid  corresponds  to  r=l  and  contains  a  3x1  transportation  problem.  Recall 
that  the  spacing  ratio  between  each  grid  in  multigrid  methods  is  2:1;  the  same 
convention  is  used  here. 

The  analogies  to  restriction,  interpolation,  and  relaxation  developed  during 
this  research  are  presented  next. 

B.  RESTRICTION 

The  restriction  method  developed  transforms  vectors  containing  fine  grid 
demands  into  vectors  containing  coarse  grid  demands,  that  is,  d^=f“d’’.  The 
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coarsest  grid  transportation  problem  is  the  3x1  problem.  This  problem  is  trivial 
to  solve  because  the  single  demand  node  receives  all  the  supply  contained  in  each 
supply  node. 

To  restrict  from  the  fine  grid  with  (2'-lf  demand  nodes  to  the  next 
coarser  grid,  with  ('2'’"”-!  demand  nodes  requires  a  method  for  aggregating 
the  demand  nodes  to  achieve  a  smaller  problem.  Let  array  dh  contain  the  fine  grid 
demand  nodes  and  array  dlh  contain  the  next  coarser  grid  demand  nodes.  The 
demand  node  d2h(i,j,k)  on  the  r"*  grid  is  made  up  of  a  weighted  average  of 
demand  node  dh(2*i,2*j,2'*k)  and  the  26  demand  nodes  surrounding  it  on  the 
(r+lf  grid.  This  is  computed  as  follows: 

d2hii,jjc)  =dh{2M,2^^^)  *i\/2)[dhi2m-h2i,2^)  ^h{2m^^-\;2*) 

*dh(2  %2  *j,2  4  -1 )  •*d/i(2  <!  +1 ,2  <^’,2  Me)  '•dft(2  ^  +1 ,2  Me) 

MiH2M,2^,2Me*\)]  •Kl/4)[dM2*i-l,2^-l,2^)  ^h{2M-l,2i;2Me-l) 
^h(2M-\;2^*\;2Me)  ^h{2M-\,2*j,2Me*l)  ^h{2M*\,2^-l,2Me) 
*dh{2M+l,2’*j,2Me-l)  *dh{2M*\,2^-*\,2Me)  *dh{2M*\,2^ilMe*\) 
+dh{2M,2*j-\,2Me-\)  ^h{2M,2i-l,2Me*\)  ^h{2M,2i*i;2Me-l) 

*iih(2  M,2  ^  +1 ,2  Me  +1 )]  +<1  /8)[dhi2  M  -1 ,2  -1 ,2  *  -1)  *dha  M  -1 ,2  -f >1 ,2  4  -1 ) 
*dh(2M-i;2^*i;2Me*l)  *dh{2M-l^^j-l^Me*l)  +dhi2M*l^i+l^Me-l) 
*dh(2M-n;2^-\;2Me+l)  *dh{2M*i;2^*i;2Me*l)  *dh(2M^i;2^-i;2Me-l)] 

The  restriction  method  is  summarized  as:  Restrict  on  fl'’ with  (2‘-l)^  demand  nodes 
to  the  coarser  grid  Q^with  (2^‘'^^-l)^  demand  nodes. 
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C  INTERPOLATION 


The  interpolation  step  transforms  vectors,  containing  the  coarse  grid  flows, 
to  vectors  containing  the  fine  grid  flows.  That  is,  x‘’=/^x^.  The  interpolation  uses 

a  local  solver  to  apportion  the  coarse  grid  flows.  Each  demand  node  on  the  coarse 
grid  represents  27  demand  nodes  on  the  fine  grid.  Thus  for  each  coarse-grid 
demand  node  we  must  apportion  the  coarse  grid  flows  among  the  27  fine  grid 
demand  nodes  in  the  least  expensive  manner.  This  is  done  by  treating  the  3 
coarse  grid  flows  to  a  given  coarse  grid  demand  node  as  supplies,  and  solving  a 
3x27  transportation  problem.  This  means  the  global  "solution"  on  is 
interpolated  to  by  solving  the  one  local  3x27  problem  for  each  demand  node 
d2h(i,j,k),  using  the  flow  Xjf2Uijjc)  supply  s,,  the  flow  as  supply  s^,  and  the 

flow  X3(2ijiM  as  supply  Sj.  The  global  solution  on  is  the  combination  of  all  the 
local  solutions.  The  3x27  problems  are  what  is  referred  to  as  the  local  transportation 
problems.  There  are  (2*''^’-l)^  local  3x27  transportation  problems  on  the  r*''  grid.  The 
local  problems  together  comprise  the  global  problem  on  the  grid  (see  Figure  11). 

As  an  example,  suppose  that  the  demand  at  node  (3,1,2)  on  a  coarse  grid  is 
for  50  units  of  a  particular  commodity,  and  that  the  coarse  grid  solution  indicates 
that  it  should  receive  40  units  from  supply  node  two,  10  units  from  supply  node 
three,  and  zero  units  from  supply  node  one.  These  flows  computed  for  demand 
node  (3,1,2)  now  become  supplies  to  the  27  associated  demand  nodes  on  the  next 
finer  grid.  Hence,  this  local  3x27  problem  has  supplies  s,=0,  $2=40,  and  $3=10. 
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Figure  11.  (a)  The  global  transportation  problem  on  a  regular  grid,  (b)  The  local 
transportation  problem. 


The  demands  for  the  local  problem  on  fl'’  must  also  be  determined.  Each 
fine  grid  demand  node  may  be  part  of  one,  two,  four,  or  eight  local  problems, 
depending  on  whether  it  lies  in  the  interior,  on  a  side,  an  edge,  or  a  comer  of  a 
local  problem.  The  total  demand  on  a  fine  grid  node  must  be  apportioned  equally 
among  all  the  local  problems  that  contain  the  node.  This  will  ensure  that  local 
demand  nodes  that  are  included  in  two  or  more  3x27  problems  do  not  receive 
more  demand  than  required  for  the  global  demand  node.  The  demand  nodes  on 
ii*’  receive  some  amount  of  demand  from  each  of  the  demand  nodes  on  Q'*. 

The  interpolation  process  can  be  summarized  as  follows: 

1.  Determine  local  demands  on  for  each  local  3x27  problem. 

2.  Solve  each  local  3x27  problem  with  the  local  solver,  using  the  flows  on  the 
coarse-grid  as  supply  values,  and  combine  the  local  solutions. 

D.  LOCAL  SOLVER 

Given  a  solution  on  0“*,  it  is  interpolated  by  solving  as  many  3x27  problems 
as  are  necessary.  The  coarsest  grid  is  a  3x1  problem,  and  only  requires  a  single 
3x27  local-solve  for  the  first  interpolation. 

Because  the  problem  is  posed  on  a  regular  grid,  there  are  thirteen  distinct 
ways  in  which  the  three  shipping  costs  for  a  demand  node  (denoted  by  i,j,  and 
it)  can  be  related.  They  are  as  follows: 
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1. 

i=;=lc 

8.  t<j<k 

2. 

i=j,  j>k 

9.  i<k<j 

3. 

H'  ;<* 

10.  j<i<k 

4. 

i=k,  j>k 

11.  j<k<i 

5. 

i=k,  j<k 

12.  k<i<j 

6. 

j=k,  j>i 

13.  k<j<i 

7. 

j=k,  j<i 

The  local  solver  exploits  these  patterns  by  choosing  an  order  in  which  to  fill  the 
demands,  using  thirteen  simple  tests.  Note  that  this  is  done  a  priori,  in  other 
words  the  list  of  ordered  choices  is  hard-coded.  This  implies  that  none  of  this 
need  be  computed;  it  is  free.  If  a  sort  routine  were  used  the  complexity  would 
grow  exponentially  with  the  size  of  the  problem;  therefore,  this  method  is  much 
more  efficient.  The  values  of  i,j,  and  Jt  relative  to  one  another  determine  the  order 
in  which  the  demand  nodes  will  be  supplied.  The  ordering  is  determined  as 
follows: 


1.  For  each  demand  node,  determine  the  cost  differential  between  the 
cheapest  (minimum(i,;,Jt))  and  the  second  cheapest  supply  node.  This 
differential  is  denoted  by  A,. 

2.  For  each  demand  node,  determine  the  cost  differential  between  the  second 
cheapest  and  the  remaining  supply  node.  The  differential  is  denoted  by  Aj. 

3.  Form  the  string  AjAj. 
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4.  Lexicographically  order  the  strings  from  largest  to  smallest.  In  the  case  of 
a  tie  the  demand  node  with  higher  overall  cost  iS  placed  first. 

The  supply  nodes  are  then  placed  in  order  of  preference  for  each  demand  node. 

The  first  choice  supply  node  will  be  the  cheapest  supply  node,  the  second  choice 

will  be  the  second  cheapest  and  the  third  choice  will  be  the  remaining  supply 

node.  The  local  solver  now  fills  each  demand  node  in  the  predetermined  order 

(the  order  given  by  the  lexicographic  ordering  of  the  strings).  For  each  demand 

node  the  method  proceeds  as  follows: 

1.  Check  the  first  choice  supply  node  for  adequate  supply.  If  there  is  more 
supply  than  demand  then  the  demand  is  filled  entirely. 

2.  If  there  is  more  demand  than  supply  and  the  supply  does  not  equal  zero, 
then  fill  the  demand  as  much  as  possible  from  the  first  choice  supply  node. 

3.  If  demand  remains,  check  the  second  choice  supply  node.  If  possible,  fill 
the  remaining  demand;  if  not,  fill  as  much  of  the  demand  as  possible  from 
the  second  choice  supply  node. 

4.  If  the  demand  is  still  not  satisfied  then  fill  the  remaining  demand  from  the 
third  choice  supply  node. 

This  local  solver  fills  each  demand  node  requirement  completely  before  repeating 
the  process  on  the  next  demand  node.  The  local  solver  provides  an  optimal 
solution  to  a  3x27  transportation  problem. 

Theorem  1,  Let  x  be  the  vector  of  flows  determined  by  the  algorithm  above  for  the 
3  X  27  problem.  Then  x  is  the  optimal  solution. 

Proof:  Let  x  be  the  flow  vector  determined  by  the  algorithm,  z  =  cx,  and  let 
z*  =  cx*  be  any  optimal  solution.  Suppose  that  x  ^  x*.  Without  loss  of  generality 
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assume  that  c,,  <  Cj,  <  Cj,.  Let  be  the  first  demand  for  which  the  flows  x  and  x* 
differ,  i.e.,  for  all  1  <  i  <  3  and  \  <  k  <  j  we  have  x*  =  x*.  As  a  consequence  of 

this  and  of  the  fact  that  the  problem  is  balanced,  for  each  1  <  i  <  3  we  have 
27 

ix^-x^  =0,  from  which  we  see  that 

ki 


17 


(6) 


3  3 

balance  implies  also  that,  for  any  choice  of  k,  ^ x*,  so  for  any  i  and  k 

i=l  i»l 

we  have 

3 

m=l 


If  we  let  A=Xj.-jc*,  then  in  particular  (7)  implies  A  .  Note  that,  since 

i  ^ 

thealgorithm  maximizes  x^,  then  A  >  0.  Now 

z*=2-tz*-z 


•z-tcx  -cx 
=z-*c{x*-x) 

3 


I-l 


;-l  27 


(8) 


By  our  choice  of  j,  the  first  sum  within  the  brackets  in  (8)  vanishes,  so  now 
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(9) 


3 

i-l 


27 


jt-y+l 

3  27 

*T.  E 

j=ijt^+i 

3  27 

"E  E 

I=U=7+1 

where  the  equality  of  (9)  and  (10)  follows  from  (7).  Rearranging,  we  have 

3  27 

2  =2-(C  -C  )A-*(C  -C3  )(x  -*3)+J)  ^ 


(10) 


(11) 


i=lJt=^+l 
27 

=2  “(c -Cj^)A -KCj^ -Cy^iXj. ~x^  +  ^  "*3*) 

fc»5+l 

27  27 

k^+\ 

-Z-(c,^-C^)A^(c^-C3;(X3^-«37 
27  27 


(12) 


(13) 


ki+l 


k^*l 


where  the  equality  of  (11)  and  (12)  follows  from  (6).  Let  5^  =  Ic,*  -  Cj^l,  and  let 
Jk-  I  ^2*  ■  ^3*  I  •  By  hypothesis,  5y  =  Cj^  -  c,^  and,  by  the  action  of  the  algorithm,  if 
k  ^ ;,  then  5^  >  5*.  For  arbitrary  values  of  Jt,  the  relationship  between  yj  and  y*  is 
not  clear,  so  let  q  =  min{-Yy,  -y*}.  Then  by  rewriting  (13)  we  have 
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27  27 

By  (6)  we  have  A=  (x*)  and  ^*“*3^)=  £  ^^3* "*3]^'  whence 
k^*[  k^+\ 


27 

2*>z-^(5^-6  )-ni  5^  ((x3j^-a:^-<X3^-x^) 
ki*\ 

=z, 

and  we  conclude  that  z  =  cx  is  optional.  QED. 

The  next  step  is  to  try  to  improve  the  answer  after  each  interpolation  sweep 
before  moving  to  the  next  finer  grid.  This  is  necessary  because,  even  though  the 
local  3x27  problems  on  the  1**'  grid  are  each  solved  to  optimality,  the  global 
problem  on  the  r""  grid  probably  won't  be  optimal.  The  reason  for  this  is  as 
follows. 

Recall  that  an  optimal  solution  to  a  transportation  problem  can  always  be 
found  containing  n+m-I  arcs.  Interpolating  from  the  coarsest  grid  (3x1)  requires 
solving  a  single  3x27  problem.  In  the  solution  to  that  problem,  no  more  than  29 
arcs  have  flow  on  them,  hence  at  most  two  demand  nodes  receive  supply  from 
more  than  one  supply  node.  As  these  flows  are  interpolated  to  the  next  finer  grid 
this  property  may  cause  difficulty.  To  interpolate,  the  local  solver  proceeds  to  fill 
the  demands  of  the  fr-I)’'  grid  local  problem  in  the  appropriate  order.  For 
example,  suppose  a  demand  node  on  the  r^grid  is  only  supplied  by  supply  node 
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one.  A  demand  node  in  the  local  problem  associated  with  this  coarse  grid  node 
whose  first  choice  supply  node  is  two  and  second  choice  supply  node  is  three  can 
only  receive  supply  from  supply  node  one.  This  is  clearly  an  optimal  solution  for 
the  local  3x27  problem  because  there  is  only  one  supply  node  containing  supply, 
but  this  is  unlikely  to  be  an  optimal  choice  in  the  global  problem.  Hence,  a 
relaxation  scheme  to  improve  the  initial  flow  allocation  before  interpolating  to  the 
next  finer  grid  is  required.  This  relaxation  scheme  is  described  next. 

E.  RELAXATION 

The  relaxation  method  checks  the  flows  computed  during  each  interpolation 
sweep  to  determine  if  each  demand  node  was  filled  by  the  cheapest  possible 
supply  node.  If  a  demand  node  is  not  supplied  by  the  cheapest  supply  node  then 
the  flow  is  diverted  to  a  temporary  demand  node,  and  the  flow  to  the  original 
demand  node  is  set  equal  to  zero.  Also,  an  equal  amount  of  flow  is  added  to  a 
temporary  supply  node,  such  that  if  the  misapportioned  supply  was  from  supply 
node  j,  the  i'*  temporary  supply  node  receives  the  same  amount  of  supply.  If  a 
demand  node  is  supplied  by  the  cheapest  supply  node  then  the  demand  in  the 
corresponding  temporary  demand  node  is  set  equal  to  zero  and  no  flow  is  added 
to  the  temporary  supply  node.  Once  all  the  demand  nodes  have  been  checked,  the 
temporary  demand  and  supply  nodes  become  part  of  a  "dummy"  problem  that 
has  the  same  structure  as  the  original  problem  but  has  only  the  supply  that  was 
misapportioned  and  has  demand  only  at  those  nodes  where  the  cheapest  supplier 
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was  not  used.  The  local  solver  is  applied  to  this  global  problem  by  solving  as 
many  local  3x27  problems  as  necessary.  The  relaxation  scheme  is  summarized  as 
follows: 


1.  If  each  demand  node  is  supplied  by  the  cheapest  available  supply  node, 
go  to  step  3;  else,  go  to  step  2  . 

2.  Divert  the  flow  to  a  temporary  demand  node  and  add  the  same  amount 
of  flow  to  a  temporary  supply  node.  Set  the  flow  to  the  original  demand 
node  equal  to  zero.  Go  to  step  4. 

3.  Set  the  corresponding  temporary  demand  node  equal  to  zero  and  do  not 
add  to  the  temporary  supply  node. 

4.  Compute  the  flows  for  the  temporary  supply  and  demand  nodes. 

5.  Add  the  new  flows  to  the  existing  solution.  Repeat  after  each  interpolation 
sweep. 

This  method  improves  the  global  flows  because  it  only  recomputes  the  flows  for 
demand  nodes  not  supplied  by  the  cheapest  supply  node  and  it  can  now  choose 
from  all  three  supply  nodes.  In  addition,  it  does  not  have  any  affect  on  the 
demand  nodes  already  supplied  by  the  cheapest  supply  node. 

The  four  operations  just  described:  restriction,  interpolation,  solve,  and 
relaxation,  are  combined  to  form  a  V-cycle  type  of  algorithm.  The  complete 
algorithm  and  the  results  of  the  tests  performed  are  presented  in  the  next  chapter. 
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V.  THE  MODIHED  V-CYCLE 


A.  DESCRIPTION  OF  THE  ALGORITHM 

In  this  chapter,  the  components  presented  in  Chapter  IV  are  combined  to 
form  a  V-cycle  algorithm.  Assume  that  there  are  t>l  grids.  The  grid  spacing  on 
the  second  coarsest  grid,  will  equal  qh  where  q=2*'^.  The  coarsest  grid 

consists  of  a  single  demand  node,  which  necessarily  receives  all  the  supply.  The 
algorithm  can  now  be  summarized  as  follows: 

1.  Restrict  the  original  3xn  transportation  problem  posed  on  O*'  to  obtain  the 
demands  on  all 

2.  Solve  the  problem  on  and  interpolate  the  solution  to  D*'. 

3.  Solve  the  problem  on  by  restricting  to  then  solving  and 
interpolating  back  to 

4.  Solve  the  problem  on  D***  by  restricting  to  D®*  then  solving  and 
interpolating  back  to  fl®*. 

5.  This  process  continues  until  fl'*’  is  reached  where  the  problem  is  solved 
and  interpolated  back  to 

6.  Compute  the  value  of  the  objective  function  on  D*’. 
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This  process  can  be  summarized  recursively  as 


O*  Solve 

Coarsen^^^^intarpoiate^Ralax 

Coaraon^—^-interpolate-Ralax 

•••  1 

O'*  Coafsen^Solve^Interpolate-Rolax 


The  net  effect  is  a  V-cycle  depicted  in  Figure  12. 


Figure  12.  The  V-Cycle  Algorithm 


This  algorithm  differs  from  the  traditional  multigrid  V-cyde  mainly  in  that 
there  is  no  analogue  to  the  residual  equation.  Thus  the  algorithm  computes  a 
solution  of  the  coarse  grid  and  interpolates  it  to  the  fine  grid,  rather  than  using 
the  coarse  grid  to  correct  the  fine  grid  solution.  The  relaxation  scheme  developed 
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here  may  provide  the  bridge  to  developing  a  FMV  cycle,  by  giving  a  mechanism 
that  can  be  used  as  a  correction  scheme.  This  will  be  discussed  in  Chapter  VI. 

This  algorithm  also  has  some  key  differences  from  aggregation  methods. 
One  difference  is  that  aggregation/disaggregation  is  performed  on  only  two 
levels;  this  algorithm  is  employed  recursively  on  several  levels.  At  first  glance  it 
may  seem  that  this  algorithm  is  increasing  the  work  required  to  solve  the 
problem,  but  this  isn't  necessarily  so.  Recall  that  on  each  level  the  problem  gets 
smaller  until  it  reaches  a  point  where  it  is  easy  and  computationally  cheap  to 
solve.  This  solution  is  then  interpolated  back  up  to  the  fine  grid/original  problem. 
This  property  is  always  exploited  in  multigrid  so  that  a  full  multigrid  V-cycle 
typically  requires  no  more  effort  than  two  fine-grid  relaxation  sweeps.  Another 
difference  is  that  aggregation  methods  only  use  a  subset  of  the  data  of  the  original 
problem  to  arrive  at  an  answer;  this  algorithm  uses  all  the  information  available. 
In  other  words,  every  demand  node  is  individually  filled  during  interpolation  and 
then  checked  for  optimality  during  relaxation.  Finally,  this  algorithm  uses  the 
flows  computed  on  the  coarser  grid  to  be  supplies  on  the  next  finer  grid,  and 
there  is  no  analogy  to  this  in  aggregation/disaggregation  methods. 

This  algorithm  is  coded  in  FORTRAN  and  each  component  is  written  as  a 
subroutine.  The  actual  code  consists  of  the  driver  code  and  four  subroutines; 
Coarsen,  Interp,  Lsolve  and  Relax.  The  data  structure  used  in  this  algorithm  is 
similar  to  that  used  in  traditional  multigrid  algorithms  and  it  is  outlined  next. 
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B.  DATA  STRUCTURE 


The  arrays  containing  the  supplies,  demands  and  flows  are  the  only 
permanent  storage  required.  All  other  calculations  are  stored  in  temporary  arrays 
which  are  overwritten  on  each  iteration.  The  supply  array  requires  three  storage 
locations  and  the  entries  remain  constant.  The  flow  array  must  be  three  times 
longer  than  the  demand  array  because  flows  from  all  three  supply  nodes  to  each 
demand  node  must  be  stored.  Demands  and  flows  are  computed  on  each  grid 
level;  therefore,  the  demands  and  flows  on  each  grid  must  be  stored. 

The  data  structure  used  stores  the  demands  and  flows,  on  each  grid 
contiguously,  in  two  arrays,  one  for  demands  and  one  for  flows.  Assume,  for 
purposes  of  this  argument,  that  the  coarsest  grid  contains  a  single  demand  node. 
Thus  a  problem  with  five  grids  contains  (2®-l)^  demand  nodes  on  the  finest  grid, 
Q**;  (2^-1)^  demand  nodes  on  (2^-1)’  demand  nodes  on  D**’;  (2^-1)^  demand 
nodes  on  Q*,  and  finally,  a  single  demand  node  on  fl’**'.  The  total  number  of 
storage  locations  required  for  the  demands  is  the  sum  of  the  storage  required  for 
the  demands  on  each  grid.  The  total  number  of  storage  locations  required  for  the 
flows  is  three  times  that  number.  For  f>l  grids,  rh  will  represent  the  grid  spacing 
on  the  coarsest  grid,  where  r=2'  *.  In  general,  the  finest  grid  requires 

(2‘-l)^+3(2'-l)^+3 

storage  locations.  Since  the  supply  storage  requirements  do  not  change  from  grid 
to  grid,  they  are  not  included  in  the  storage  calculations.  For  ease  of  calculation 
let  N=2‘  be  the  upper  bound,  since  2‘-l<2'.  Thus  D'"  requires  slightly  fewer  than 
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4N^  storage  locations  for  demands,  flows,  and  supplies.  Moving  from  any  grid  to 
the  next  coarser  grid  reduces  the  number  of  storage  locations  by  a  factor  of  2^=8. 
Adding  the  storage  requirements  on  each  grid  and  using  the  sum  of  the  resulting 
geometric  series  yields 

Sforflge=4N'(l+2-"+2^+...+2'^)<4NV(l-2-')=(32/7)N' 
which  is  an  upper  bound  on  permanent  storage  locations. 

C  COMPUTATIONAL  COMPLEXITY 

The  driver  code  and  the  solve  subroutine  both  have  constant  complexity, 
denoted  by  0(1).  For  f>l  grids,  rh  represents  the  grid  spacing  on  the  coarsest  grid, 
where  r=2'  ’.  The  Coarsen  subroutine  has  a  constant  number  of  calculations  per 
grid  point,  denoted  by  /.  These  calculations  are  performed  (2‘'  '’-l)^  times  on  the 
finest  grid,  £1^.  These  calculations  are  then  performed  (2"’^-l)^on  O'",  and  so  on 
until  the  grid  preceding  the  coarsest  grid  is  reached,  where  these  operations 
are  performed  (2^-1)^  times.  For  ease  of  calculation,  let  M=2*  ’  be  the  upper  bound, 
since  Therefore,  £2'’  requires  /M’  operations.  Moving  from  any  grid  to 

the  next  coarser  grid  reduces  the  number  of  operations  by  a  factor  or  2^=8.  This 
gives  a  total  operation  count  =  /M^(l+2'’+2'*+...+2'^)</M’/(l-2'^)  =  (8/7)/M^ 

A  similar  argiuuent  can  be  given  for  the  Interp  subroutine.  Let  K  denote  the 
number  of  constant  calculations  per  grid  point  in  Interp;  therefore  the  total 
number  of  operations  in  Interp  is  given  by  Tf={8/7)KM^.  For  the  Relax  subroutine 
the  same  argument  is  used  again,  but  in  this  case  N=2'  is  the  upper  bound.  This 
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is  because  the  relaxation  scheme  is  performed  at  most,  (2‘)^  times,  on  any  grid. 
Coarsen  and  Interp  are  performed  at  most  (2'  ')^  times  on  any  grid.  Let  L  denote 
the  number  of  constant  calculations  required  in  Relax.  Then  the  total  number  of 
calculations  required  in  Relax  is  TR=(8/7)LN^  Considering  the  worst-case  scenario, 
and  the  fact  that  N>M;  the  computational  complexity  of  the  entire  algorithm  can 
be  expressed  as  S(8/7)N^  where  S=/+jK+L.  Therefore,  this  algorithm  has  O(N^) 
complexity. 

D.  PERFORMANCE 

A  separate  program  was  written  to  generate  random  supply  and  demand 
node  quantities,  which  were  used  to  generate  test  problems.  The  supplies, 
demands  and  costs  of  the  test  problems  were  all  integer  valued.  Approximately 
twenty  sample  transportation  problems  with  random  demand  and  supply 
quantities  were  solved  using  3,4  and  5  grid  levels.  The  average  results  on  each 
grid  are  depicted  in  Table  II.  Table  II  shows  that  this  algorithm  comes  consistently 
within  4%  of  optimality. 

These  results  indicate  that  multigrid  techniques  can  be  successfully  used  to 
solve  large  scale  long  transportation  problems  posed  in  cost-space  on  regular 
grids.  The  solutions  obtained  are  fairly  close  to  optimality.  While  these  results  are 
very  encouraging,  there  remain  some  difficult  aspects  of  the  algorithm  to  be 
addressed. 
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First,  the  algorithm  produces  real-valued  demands  and  flows  through  the 
restriction  and  interpolation  schemes.  However,  many  real  world  problems  of  this 
type  require  integer  output. 

A  second  problem  is  that  some  flow  is  lost  due  to  round-off  error.  Therefore, 
total  supply  does  not  equal  total  flow  output.  The  flow  loss  increases  with  the 
size  of  the  problem.  For  example,  the  3x343  problem  loses  no  measurable  amount 
of  flow.  The  3x3375  problem  loses  .0003%  of  the  total  flow.  And  finally,  the 
3x29,791  problem  loses  .0026%  of  the  total  flow.  This  problem  may  disappear  with 
improved  restriction  and  interpolation  schemes.  These  improved  schemes  must 
be  able  to  preserve  the  integer  values  of  the  input  to  provide  integer  output. 

A  third  problem  stems  from  the  fact  that  there  are  more  than  m+n-1  arcs 
having  flow  on  them  in  the  final  solution,  which  is  an  indicator  that  the  solution 
may  not  be  optimal.  For  example,  the  3x343  problem  should  contain,  at  most,  345 
arcs  with  nonzero  flow  in  the  final  solution,  but  there  are  actually  423  such  arcs 
in  the  final  solution.  This  implies  that  the  underlying  undirected  graph  of  the 
solution  to  the  transportation  problem  contains  cycles.  The  interpolation  scheme 
used  in  this  algorithm  is  responsible  for  introducing  these  cycles.  The  cycles  are 
a  result  of  the  coarse  grid  flows  supplying  only  a  portion  of  the  demand  for  each 
demand  node  on  the  next  finer  grid.  Thus  each  demand  node  on  the  next  finer 
grid  could  ultimately  be  supplied  by  up  to  three  different  supply  nodes.  A  search 
method  to  detect  and  eliminate  cycles  would  be  computationally  expensive 
(Cavanaugh,  1992).  However,  revisions  to  the  relaxation  scheme  used  in  this 
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algorithm  could  prove  beneficial.  Some  suggestions  are  provided  m  the  next 
chapter,  as  are  conclusions,  and  recommendations  for  future  research. 


TABLE  II.  V-CYCLE-TYPE  ALGORITHM 


Number  of  Grids 

3xn 

Algorithm  Time  to  Solve 

%  Over  Optimal 

2 

27 

0.001s 

0 

3 

343 

0.63s 

3.8 

4 

3,375 

8.0s 

3.4 

5 

29,791 

78.41s 

4.0 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

This  thesis  shows  that  multigrid  methods  can  be  successfully  applied  to 
solution  of  large  scale  transportation  problems  when  they  are  posed  on  a  regular 
grid  in  cost  space.  The  method  developed  here  is  not  intended  to  replace  any 
existing  solution  method.  Instead  the  original  goal  of  this  thesis  is  to  determine 
the  applicability  of  the  multilevel  approach  to  solving  the  long  transportation 
problem  posed  on  a  regular  grid.  The  initial  results  are  very  encouraging,  but 
more  research  is  needed. 

A  significant  contribution  of  this  thesis  is  the  identification  of  a  relaxation 
scheme.  The  relaxation  method  developed  is  able  to  correct  flows  locally  that  are 
not  supplied  by  the  cheapest  supply  node.  These  corrections  are  made  locally  on 
the  global  grid  on  each  grid  level  after  the  interpolation  scheme.  A  method  similar 
to  this  has  not  yet  been  developed  for  solving  transportation  problems  using 
multigrid  methods.  In  addition,  this  method  could  lead  to  a  cycle  detection 
scheme  so  that  the  solution  would  contain  at  most  m+n-1  arcs  in  the  final 
solution. 

The  relaxation  scheme  may  also  be  the  key  to  developing  a  full  multigrid  V- 
cycle.  Call  the  misdirected  flow,  that  is  placed  in  the  temporary  array,  the  error. 
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Instead  of  solving  several  3x27  local  problems  to  correct  the  flow,  solve  one  global 
problem.  The  global  problem  is  of  the  form  Ae=r.  The  question  now'  arises:  How 
does  one  solve  Ae=r?  The  answer  is  to  use  the  V-Cycle  algorithm  de\'eloped  in 
this  thesis.  Recall  that  u  can  be  replaced  by  e  and  f  can  be  replaced  bv  r  in  tne 
equation  Au=f  to  arrive  at  an  equation  of  the  same  form,  Ae=r.  Suppose  there 
exists  an  initial  guess  on  the  coarse  grid,  interpolate  the  fir  w  to  the  next  level  and 
find  the  error  using  the  relaxation  scheme.  Repeat  until  the  fine  grid  answer  is 
reached.  This  algorithm  results  in  a  Full  Multigrid  V-cycle  algorithm  depicted  in 
Figure  9(b). 

The  problem  solved  in  this  thesis  is  deliberately  arhficial.  Further 
advancements  in  this  area  of  research  require  a  soluhon  method  that  is  applicable 
to  more  general  problems.  Specifically,  a  solver  that  can  handle  more  than  three 
supply  nodes  must  be  developed. 

B.  AREAS  FOR  FUTURE  RESEARCH 

The  results  of  this  thesis  indicate  that  further  research  in  this  area  can  be 
productive.  There  are  several  avenues  for  future  research. 

This  algorithm  solves  a  very  specific  problem  that  is  posed  on  a  three- 
dimensional  grid.  Generalization  to  an  m-dimensional  grid  would  logically  follow. 
Research  into  developing  a  method  to  map  an  arbitrary  transportation  problem 
to  a  regular  grid  and  then  map  the  solution  back  to  the  original  problem  domain 
is  needed. 
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Another  area  of  research  that  is  needed  for  the  algorithm  developed  is  to 
find  a  different  method  of  restriction  that  will  not  introduce  real-valued  demands 
or  flows  into  the  problem.  One  such  method  could  be  analogous  to  the  finite 
volume  element  discretization  used  in  modem  multigrid  methods  (McCormick, 
1989).  A  new  interpolation  method  that  also  preserves  integer  flows  is  necessary. 
These  are  the  two  most  urgent  areas  of  follow-on  research  aimed  at  correcting  the 
problem  of  lost  flow  in  the  output.  As  stated  earlier,  most  users  require  integer 
flows  and  costs. 

A  third  area  of  research  is  to  devise  a  more  efficient  relaxation  scheme,  or 
perhaps  to  develop  a  new  relaxation,  based  on  detecting  and  eliminating  cycles 
(Cavanaugh,  1992).  The  scheme  developed  here  already  contains  the  necessary 
information,  such  as  the  flow  for  each  demand  node  and  the  supplier.  Ar  2fficient 
means  of  using  this  data  to  find  cycles  must  be  determined.  A  test  similar  to  the 
existing  one  that  checks  for  optimality  could  also  check  for  multiple  supply  nodes 
to  a  single  demand  node  and  then  redirect  the  flow. 

Lastly  and  most  importantly,  this  thesis  provides  a  direction  to  begin  to 
implement  a  true  Full  Multigrid  V-cycle  algorithm  which  has  not  been  done 
before. 
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APPENDIX 


C  THIS  CODE  WAS  WRITTEN  BV  ANNETTE  P.  CORNETT 

C  1606  KICKAPOO  RD. 

C  PUEBLO,  CO  81001 

C 

C  THIS  PROGRAM  WILL  SOLVE  A  LONG  TRANSPORTATION  PROBLEM  ON  A  REGULAR 

C  GRID  IN  COST  SPACE. THE  LARGEST  SIZE  PROBLEM  SOLVED  BY  THIS  PROGRAM 

C  IS  A  3X29,071  TRANSPORTATION  PROBLEM.  THE  PARAMETER,  NGRIDS  MUST 

C  BE  CHANGED  BEFORE  THE  PROGRAM  IS  RUN.  THE  ARRAY  FLOW  CONTAINS  THE 

C  FLOWS  ON  EACH  GRID.  IT  IS  INITIALIZED  TO  ZERO  IN  THE  BEGINNING. 

C  THE  ARRAY  DEMANDS  CONTAINS  THE  QUANTITY  OF  DEMAND  FOR  EACH  DEMAND 

C  NODE  ON  EACH  GRID.  THE  ARRAY  SUPPLY  CONTAINS  THE  INITIAL  SUPPLY 

C  FOR  THE  FINE  GRID  PROBLEM.  THIS  PROGRAM  CONTAINS  THREE 

C  SUBROUTINES,  THEY  ARE:  COARSEN,  INTERP,  LSOLVE  AND  RELAX.  A  BRIEF 

C  EXPLANATION  FOR  EACH  PRECEDES  THE  ROUTINE. 

PROGRAM  TRANS 11 
INTEGER  DIM 

PARAMETER  (NGRIDS  =  3, DIM  =  3) 

REAL  FLOW(90000) , DEMANDS (30000) ,SUPPLY(DIM) , XHTEMP (DIM, 63 , 63 , 

+63 ) , DHTEMP (63,63,63), STEMP (DIM) 

INTEGER  NROOT(NGRIDS) ,IDX(NGRIDS) ,ED(NGRIDS) , NPTS (NGRIDS ) ,L, COUNT 
COUNT  =0 
ITEMP  =  1 
JTEMP  =  1 
KTEMP  =  1 
DO  1  1=) ,3 
1  STEMP (I) =0 

C  SET  UP  INDEX  SCHEME  FOR  ARRAYS  DEMANDS  AND  FLOWS 

DO  5  I  =  1, NGRIDS 
NROOT(I)  =  2**I-1 
NPTS(I)  =  (2**I-1) **DIM 
IDX(I)  =  1 
ED ( I )  =  1 
5  CONTINUE 

DO  10  J  =  2, NGRIDS 

IDX(J)  =  IDX(J-l)  +  NPTS(J-l) 

ED(J)  =  IDX(J)  +  NPTS(J)  -1 
10  CONTINUE 

K  =  IDX (NGRIDS)  +  NPTS (NGRIDS)  -  1 
DO  15  I  =  1,K 

DEMANDS ( I )  =0 
15  CONTINUE 

DO  20  I  =  1,3*K 
FLOW (I)  =  0 
20  CONTINUE 

OPEN(UNIT  =  6, FILE  =  'FINE  DATA  D'. STATUS  =  'OLD') 

READ (6, 22) (SUPPLY (I) , 1=1, DIM) 

22  FORMAT (F12.0) 

READ(6,22) (DEMANDS(I) , I=IDX (NGRIDS) ,ED(NGRIDS) ) 

DO  25  I  =  NGRIDS, 3,-1 
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25  CALL  COARSEN (DEMANDS ( IDX ( I ) ) , DEMANDS ( IDX ( I-l ) ) , NROOT ( I ) , NROOT ( I -1 ) 
+  ) 

1=2 

CALL  LSOLVE(FLOW(3*IDX(I) -2) , DEMANDS ( I DX ( I ) ) .SUPPLY, DIM, 

+ ITEMP , JTEMP , KTEMP ) 

DO  30  I  =  2,NGRIDS-1 
DO  31  K  =  l.NROOTd) 

DO  31  J  -  l.NROOTd) 

DO  31  II  =  l.NROOTd) 

DO  31  L  =  l.DIM 

XHTEMP(L.2*II.2*J-2*K)  =  0 

IF( (K.LE.3) .AND. (J.LE.3) .AND. dI.LE.3) .AND. (L.EQ.l) 

+ ) THEN 

DHTEMP (11 , J,K)  =  0 
ENDIF 

31  CONTINUE 

CALL  INTERP (DEMANDS (IDX (I +1) ) , DHTEMP , XHTEMP , STEMP , 
+FLOW(3*IDX(I) -2) .FLOW (3* IDX (1+1) -2) , NROOT (I ) .NROOT (1+1) .DIM) 

30  CALL  RELAX(STEMP,FLOW(3*IDXd  +  l) -2)  .NROOT(I)  ,NROOT(I  +  l)  .DIM) 

END 

C  THE  SUBROUTINE  COARSEN  WILL  RESTRICT  THE  FINE  GRID  DEMAND  NODES  TO 

C  THE  NEXT  COARSER  GRID  BY  USING  A  WEIGHTED  AVERAGE.  THE  ARRAY 

C  DH  CONTAINS  THE  DEMAND  FOR  EACH  FINE  GRID  DEMAND  NODE.  THE 

C  ARRAY  D2H  CONTAINS  THE  DEMAND  FOR  EACH  FINE  GRID  DEMAND  NODE, 

C  R  IS  THE  ROOT  OF  THE  FINE  GRID,  RC  IS  THE  ROOT  OF  THE 

C  COARSE  GRID. 

SUBROUTINE  COARSEN (DH , D2H , R, RC) 

INTEGER  R.RC.SUM 

REAL  IM, JM.KM, IP, JP.KP, IMJM, IMJP, IMKM, IMKP, JMKM. JMKP. IPJP, IPJM. 
+IPKP, IPKM, JPKP, JPKM, IMJMKM, IMJPKM. IMJMKP. IMJPKP, IPJPKP, IPJMKM, 
+IPJPKM, IPJMKP,DH(R,R,R) , D2H (RC. RC, RC) 

SUM  =  0 

C  COMPUTE  THE  WEIGHTED  AVERAGE 


DO  10  K  =  l.RC 

DO  10  J  =  1,  RC 
DO  10  I  =  l.RC 


D2H(I 

,  J 

,K) 

IM  = 

.5 

IP  = 

.5 

JM  = 

.5 

JP  = 

.5 

KM  = 

.5 

KP  = 

,5 

IMJM 

= 

.25 

IMJP 

= 

.25 

IPJM 

= 

.25 

IPJP 

= 

.25 

IMKM 

= 

.25 

IMKP 

= 

.25 

IPKM 

= 

.25 

IPKP 

= 

.25 

JMKM 

= 

.25 

JMKP 

= 

.25 

JPKM 

= 

.25 

D2Hd,J,K) 


DH(2*I,2*J,2*K) 
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JPKP  =  .25 
IMJMKM  =  .125 
IMJPKM  =  .125 
IMJMKP  =  .125 
IMJPKP  =  .125 
IPJMKM  =  .125 
IPJMKP  =  .125 
IPJPKM  =  .125 
IPJPKP  =  .125 

C  THE  IF  STATEMENTS  CHECK  FOR  EDGE  AND  CORNER  NODES 

IF  (I  .EQ.  1)  THEN 
IM  =  1 
IMJM  =  . 5 
IMJP  =  .5 
IMKM  =  .5 
IMKP  =  .5 
IMJPKM  =  .25 
IMJMKP  =  .25 
IMJPKP  =  .25 
IMJMKM  =  .25 

ELSEIFd  .EQ.RC)  THEN 
IP  =  1 
IPJM  =  .5 
IPJP  =  .5 
IPKM  =  .5 
IPKP  =  .5 
IPJPKM  =  .25 
IPJMKP  =  .25 
IPJPKP  =  .25 
IPJMKM  =  .25 
ENDIF 

IF  (J  .EQ.  1)  THEN 
JM  =  1 
IMJM  =  .5 
IPJM  =  .5 
JMKM  =  .5 
JMKP  =  .5 
IMJMKM  =  .25 
IPJMKP  =  .25 
IMJMKP  =  .25 
IPJMKM  =  .25 
ELSEIF  (J  .EQ.RC)  THEN 
JP  =  1 
IMJP  =  .5 
IPJP  =  .5 
JPKM  =  .5 
JPKP  =  .5 
IMJPKM  =  .25 
IMJPKP  =  .25 
IPJPKP  =  .25 
IPJPKM  =  .25 
ENDIF 

IF  (K  .EQ.  1)  THEN 
KM  =  1 
IMKM  =  .5 
IPKM  =  .5 
JMKM  =  .5 
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JPKM  = 

.5 

IMJPKM 

=  .25 

IMJMKM 

=  .25 

IPJPKM 

=  .25 

IPJMKM 

=  .25 

ELSEIF  (K 

.EQ.RC 

KP  =  1 

IMKP  = 

.5 

IPKP  = 

.5 

JMKP  = 

.5 

JPKP  = 

.5 

IMJPKP 

=  .25 

IMJMKP 

=  .25 

IPJPKP 

=  .25 

IPJMKP 

=  .25 

ENDIF 


THEN 


D2H(I,J,K)  =  D2H{I,J.K)  +  IM*DH ( 2 • I - 1 , 2 *J , 2*K ) 
+  IP*DH(2*l4-l,2*J.2*K)  t  JM*DH(2*1,2*J-1,2*K) 

+  JP*DH(2*I,2*J+1.2»K)+  KM*DH(2*I,2*J,2*K-1) 

+  KP*DH(2*I,2*J,2*K+1) 

IF  ((I  .EQ.l)  .AND.  (J  .EQ.  1))  THEN 
IMJM  =  1 
IMJMKP  =  .  5 
IMJMKM  =  . 5 

ELSEIF  ((I  .EQ.l)  .and.  (J  .EQ.RC))  THEN 
IMJP  =  1 
IMJPKP  =  .5 
IMJPKM  =  .5 
ENDIF 

IF  ((I  .EQ.RC)  .AND.  (J  .EQ.  1))  THEN 
IPJM  =  1 
IPJMKP  =  .5 
IPJMKM  =  .5 

ELSEIF  ((I  .EQ.RC)  .AND.  (J  .EQ.RC))  THEN 
IPJP  =  1 
IPJPKP  =  .5 
IPJPKM  =  .5 
ENDIF 

IF  ((I  .EQ.l)  .AND.  (K  .EQ.  1))  THEN 
IMKM  =  1 
IMJMKM  =  . 5 
IMJPKM  =  .5 

ELSEIF  ((I  .EQ.l)  .and.  (K  .EQ.RC))  THEN 
IMKP  =  1 
IMJMKP  =  .5 
IMJPKP  =  .5 
ENDIF 

IF  ((I  .EQ.RC)  .AND.  (K  .EQ.  1))  THEN 
IPKM  =  1 
IPJMKM  =  .5 
IPJPKM  =  .5 

ELSEIF  ((I  .EQ.RC)  .AND.  (K  .EQ.RC))  THEN 
IPKP  =  1 
IPJMKP  =  .5 
IPJPKP  =  .5 
ENDIF 

IF  ( (J  .EQ.  1)  .AND.  (K  .EQ.  1))  THEN 
JMKM  =  1 
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n  n  n  n 


IMJMKM  =  . 5 
I P JMKM  =  . 5 

ELSEIF  ( (J  .EQ.  1)  .AND.  (K  .EQ.RC))  THEN 
JMKP  =  1 
IMJMKP  =  .5 
IPJMKP  =  .5 

END  IF 

IF  ( (J  .EQ.RC)  .AND.  (K  . EQ .  1))  THEN 
JPKM  =  1 
IMJPKM  =  .5 
IPJPKM  =  .5 

ELSEIF  ( (J  .EQ.RC)  .AND.  (K  .EQ.RC))  THEN 
JPKP  =  1 
IMJPKP  =  .5 
IPJPKP  =  .5 

ENDIF 

D2H(I,J,K)  =  D2H(I,J,K)  +  IMJM«DH (2* I - 1 , 2 *J- 1 , 2*K) 

+  +  IPJM*DH{2*I  +  1,2*J-1,2«K)  +  IMJP*DH ( 2  *  I - 1 , 2* J+ 1 , 2*K) 

4-  +  IPJP*DH(2*I  +  1,2*J  +  1.2*K)  +  IMKM*DH(2*I-1,2*J,2*K-1) 

f  4-  IMKP*DH(2*I-1.2*J.2«K  +  1)  +  I  PKM'DH  ( 2  •  I  + 1 , 2  *  J  ,  2 'K- 1 ) 

♦  V  IPKP'DH (2»I  +  1 , 2*J, 2»K  +  1 )  ♦  JMKM*DH ( 2 *  I , 2 * J -  1 , 2 *K - 1 ) 

4-  +  JMKP*DH(2*I,2*J-1,2*K+1)  +  JPKM«DH  (2*1 , 2*J  +  1 , 2*K-1 ) 

+  4  JPKP*DH(2*I,2*J+1,2*K+1) 


IF  ( (I.EQ.l) .AND. (J.EQ.l) .AND. (K.EQ.l) )  THEN 
IMJMKM  =  1 

ELSEIF  ( (I.EQ.l) .AND, (J.EQ.l) .AND. (K. EQ.RC) )  THEN 
IMJMKP  =  1 

ELSEIF  ( (I.EQ.l) .AND. (J. EQ.RC) .AND. (K.EQ.l) )  THEN 
IMJPKM  =  1 

ELSEIF  ( (I ,EQ.l) .AND. (J. EQ.RC) .AND. (K. EQ.RC) )  THEN 
IMJPKP  =  1 

ELSEIF  ( (I. EQ.RC) .AND. (J.EQ.l) .AND. (K.EQ.l) )  THEN 
IPJMKM  =  1 

ELSEIF  ( (I. EQ.RC) .AND. (J. EQ.RC) .AND. (K.EQ.l) )  THEN 
IPJPKM  =  1 

ELSEIF  ( (I. EQ.RC) .AND. (J. EQ.RC) .AND. (K, EQ.RC) )  THEN 
IPJPKP  =  1 

ELSEIF  ( (I. EQ.RC) .AND. (J.EQ.l) .AND. (K. EQ.RC) )  THEN 
IPJMKP  =  1 

ENDIF 

D2H(I,J,K)  =  D2H(I,J,K)  4-  IMJMKM*DH  (2*1-1 , 2*J-1 , 2*K-1 ) 

4-  4-IPJMKM*DH(2*l4-l,2*J-l,2*K-l)  4-IPJPKM*DH(2*l4.1,2*J4-l,2*K-l) 

4-  4-IPJMKP*DH(2*I  +  1,2*J-1,2*K+1)  4-IMJPKM*DH(2*I-l,2*J4.1,2*K-l) 

+  •♦•IMJMKP*DH(2*I-1.2*J-1.2*K+1)  +IMJPKP*DH(2*I-1,2*J4-1.2*K+1) 

4-  4-IPJPKP*DH(2*I+l,2*J+l,2*K+l) 

10  CONTINUE 


RETURN 

END 

THE  SUBROUTINE  LSOLVE  SOLVES  THE  LOCAL  3X27  TRANSPORTATION  PROBLEM 
THE  ARRAY  X2H  CONTAINS  THE  FLOWS  FOR  EACH  COARSE  GRID  DEMAND  NODE. 
THE  ARRAY  D2H  CONTAINS  THE  DEMANDS  FOR  EACH  LOCAL  TRANSPORTATION 
PROBLEM  DEMAND  NODE.  THE  DEMAND  NODES  ARE  FILLED  IN  LEXICOGRAPHIC 
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C  DEPENDING  ON  THE  RELATIONSHIP  BETWEEN  THE  COSTS  FROM  EACH  SUPPLY 

C  NODE.  THE  ARRAY  D  CONTAINS  THE  FIRST  CHOICE  SUPPLY  NODE  FOR  EACH 

C  DEMAND  NODE.  THE  ARRAY  DP  CONTAINS  THE  SECOND  CHOICE  SUPPLY  NODE 

C  FOR  EACH  DEMAND  NODE.  THE  ARRAY  DPP  CONTAINS  THE  THIRD  CHOICE 

C  SUPPLY  NODE  FOR  EACH  DEMAND  NODE. 

SUBROUTINE  LSOLVE  ( X2H , D2H, S , DIM, I I , JJ , KK ) 

INTEGER  DIM, P,D(27) ,DP(27) , DPP (27) 

REAL  D2H(3, 3,3) ,  X2H ( DIM , DIM, DIM , DIM) , S (DIM) ,  X ( 3 , 27 ) , DEMAND ( 27 ) , 
+DTEMP (27 ) 

•  INITIALIZE  ARRAYS 

DO  1  I  =  1,27 
DTEMP ( I )  =0 
D(I)  =  0 
DP (I)  =  0 
DPP (I)  =  0 
DEMAND ( I )  =0 
1  CONTINUE 

P  =  1 

DO  5  L  =  1,DIM 
DO  5  K  =  1,3 

DO  5  J  =  1,3 
DO  5  I  =  1,3 

X2H(L,I,J,K)  =  0 
IF  (L  .EQ.  DIM)  THEN 

DTEMP (P)  =  D2H(I,J,K) 

P  =  P  +  1 
ENDIF 

5  CONTINUE 

DO  6  I  =  l.DIM 
DO  6  J  =  1,27 
X(I,J)  =  0 
6  CONTINUE 


REORDER  DEMAND  NODE  ARRAY 


IF( (II .EQ, JJ) .AND. (JJ.EQ.KK) )  THEN 


DEMAND ( 1 ) 
DEMAND ( 2 ) 
DEMAND (3) 
DEMAND (4) 
DEMAND (5) 
DEMAND (6) 
DEMAND (7) 
DEMAND(8) 
DEMAND (9) 
DEMAND (10) 
DEMAND (11) 
DEMAND (12) 
DEMAND (13) 
DEMAND (14) 
DEMAND(15) 
DEMAND (16) 
DEMAND (17) 
DEMAND (18) 
DEMAND (19) 
DEMAND(20) 


=  DTEMP (9) 

=  DTEMP(21) 

=  DTEMP (25) 

=  DTEMP (6) 

=  DTEMP (8) 

=  DTEMP (12) 

=  DTEMP (16) 

=  DTEMP(20) 

=  DTEMP (22) 

=  DTEMP (18) 
=  DTEMP (24) 
=  DTEMP(26) 
=  DTEMP (5) 

=  DTEMP (11) 
=  DTEMP (13) 
=  DTEMP (3) 

=  DTEMP (7) 

=  DTEMP (19) 
=  DTEMP (15) 
=  DTEMP (17) 
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DEMAND(21)  =  DTEMP(23) 

DEMAND(22)  =  DTEMP(2) 

DEMAND (23)  =  DTEMP(4) 

DEMAND (24)  =  DTEMP(IO) 

DEMAND (25)  =  DTEMP(l) 

DEMAND (26)  =  DTEMP(14) 

DEMAND(27)  =  DTEMP(27) 

D(l)  =  3 

DP(1)  ^  1 

DPP(l)  =  2 

D(2)  =  2 

DP(2)  =  1 

DPP(2)  =  3 

D(3)  =  1 

DP(3)  =  3 

DPP(3)  =  2 

D(4)  =  3 

DP(4)  =  2 

DPP (4)  =  1 

D(5)  =  3 

DP (5)  =  1 

DPP (5)  =  2 

D(6)  =  2 

DP (6)  =3 

DPP (6)  =  1 

D(7)  =  1 

DP (7)  =3 

DPP(7)  =  2 

D(8)  =  2 

DP(8)  =  1 

DPP (8)  =  3 

D(9)  =  1 

DP(9)  =  2 

DPP(9)  =  3 

D(10)  =  3 

DP(IO)  =  2 

DPP(IO)  =  1 

D(ll)  =  2 

DP (11)  =  1 

DPP(ll)  =  3 

D(12)  =  1 

DP(12)  =  3 

DPP(12)  =  2 

D(13)  =  3 

DP (13)  =2 

DPP (13)  =  1 

D(14)  =  2 

DP (14)  =  1 

DPP(14)  =  3 

D(15)  =  1 

DP (15)  =  2 

DPP(15)  =  3 

D(16)  =  3 

DP(16)  =  2 

DPP (16)  =  1 

D(17)  =  3 

DP (17)  =  1 

DPP (17)  =  2 

D(18)  =  2 

DP(18)  =  1 
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3 


DPP(18)  = 
D(19)  =2 
DP(19)  =  3 
DPP(19)  =  1 
D(20}  =  1 
DP{20)  =  2 

DPP(20)  =  3 

D(21)  =2 
DP(21)  =  1 
DPP (21)  =  3 
D(22)  =  3 
DP(22)  =  2 
DPP(22)  =  1 
D(23)  =  1 
DP(23)  =  3 
DPP(23)  =  2 
D(24)  =  2 
DP(24)  =  1 
DPP(24)  =  3 
D(25)  =  1 
DP(25)  =  2 
DPP(25)  =  3 
D(26)  =  2 
DP(26)  =  3 
DPP(26)  =  1 
D(27)  =  3 
DP(27)  =  1 
DPP(27)  =  2 


ELSEIF( (Il.EQ.JJ) .AND. (KK.GT.JJ) ) 


DEMAND (1) 
DEMAND (2) 
DEMAND ( 3 ) 
DEMAND(4) 
DEMAND (5) 
DEMAND (6) 
DEMAND (7) 
DEMAND ( 8 ) 
DEMAND ( 9 ) 
DEMAND (10) 
DEMAND (11) 
DEMAND (12) 
DEMAND (13) 
DEMAND (14) 
DEMAND (15) 
DEMAND (16) 
DEMAND (17) 
DEMAND(18) 
DEMAND (19) 
DEMAND (20) 
DEMAND(21) 
DEMAND(22) 
DEMAND(23) 
DEMAND(24) 
DEMAND(25) 
DEMAND(26) 
DEMAND (27) 
DO  10  I  = 
D(I) 


DTEMP(21) 
DTEMP(25) 
DTEMP(12) 
DTEMP(16) 
DTEMP(3) 
DTEMP(7) 
DTEMP(20) 
DTEMP(22) 
DTEMP(24) 
DTEMP (26) 
DTEMP (11) 
DTEMP (13) 
DTEMP (15) 
DTEMP (17) 
DTEMP (2) 
DTEMP ( 4 ) 
DTEMP ( 6 ) 
DTEMP (8) 
DTEMP (19) 
DTEMP (23) 
DTEMP(IO) 
DTEMP (27) 
DTEMP (14) 
DTEMP (1) 
DTEMP (18) 
DTEMP ( 5 ) 
DTEMP (9) 
21 .2 
2 


THEN 
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DP (I)  =  1 
DPP (I)  =  3 

10  CONTINUE 

DO  11  I  =  2.26,2 
D(I)  =  1 
DP ( I )  =2 
DPP(I)  =  3 

11  CONTINUE 

ELSEIFC (lI.EQ.JJ) .AND. (JJ.GT.KK) )  THEN 

DEMAND  (1)  =  DTEMPO) 

DEMAND(2)  =  DTEMP(6) 

DEMAND (3)  =  DTEMP(8) 

DEMAND(4)  =  DTEMP(18) 

DEMAND ( 5 )  =  DTEMP ( 5 ) 

DEMAND(6)  =  DTEMP(3) 

DEMAND (7)  =  DTEMP (7) 

DEMAND{8)  =  DTEMP(15) 

DEMAND (9)  =  DTEMP (17) 

DEMAND (10)  =  DTEMP (2) 

DEMAND (11)  =  DTEMP (4) 

DEMAND (12)  =  DTEMP (27) 

DEMAND (13)  =  DTEMP (14) 

DEMAND (14)  =  DTEMP (1) 

DEMAND (15)  =  DTEMP (12) 

DEMAND(16)  =  DTEMP(16) 

DEMAND (17)  =  DTEMP (24) 

DEMAND(18)  =  DTEMP(26) 

DEMAND(19)  =  DTEMP(ll) 

DEMAND (20)  =  DTEMP (13) 

DEMAND (21)  =  DTEMP (23) 

DEMAND (22)  =  DTEMP (10) 

DEMAND(23)  =  DTEMP(21) 

DEMAND (24)  =  DTEMP (25) 

DEMAND (25)  =  DTEMP (20) 

DEMAND(26)  =  DTEMP(22) 

DEMAND (27)  =  DTEMP (19) 

DO  12  I  =  1,13,2 
D(I)  =  3 
DP (I)  =  1 
DPP (I)  =  2 

12  CONTINUE 

DO  13  I  =  2,14,2 
D(I)  =  3 
DP(I)  =  2 
DPP (I)  =  1 

13  CONTINUE 

DO  112  I  =  15,27,2 
D(I)  =  3 
DP (I)  =  2 
DPP (I)  =  1 

112  CONTINUE 

DO  113  I  =  16,26,2 
D(I)  =  3 
DP (I)  =  1 
DPP(I)  =  2 

113  CONTINUE 

ELSEIF( (II.EQ.KK) .AND. (JJ.GT.II) )  THEN 
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DEMAND  ( 1 )  =  DTEMP  ( 9  ) 
DEMAND(2)  =  DTEMP(25) 
DEMAND  ( 3  )  =  DTEMP  ( 6  ) 
DEMAND(4)  =  DTEMP(22) 
DEMAND (5)  =  DTEMP (3) 

DEMAND(6)  =  DTEMP (19) 
DEMAND  (7)  =  DTEMP  (8) 

DEMAND  (8)  =  DTEMP  (16) 
DEMAND (9)  =  DTEMP (18) 
DEMAND(IO)  =  DTEMP(26) 
DEMAND (11)  =  DTEMP ( 5 ) 
DEMAND(12)  =  DTEMP(13) 
DEMAND (13)  =  DTEMP (15) 
DEMAND (14)  =  DTEMP (23) 

DEMAND (15)  =  DTEMP ( 2 ) 
DEMAND  (16)  =  DTEMP  (10) 

DEMAND (17)  =  DTEMP (12) 
DEMAND (18)  =  DTEMP (20) 
DEMAND (19)  =  DTEMP (7) 
DEMAND(20)  =  DTEMP(17) 
DEMAND(21)  =  DTEMP(4) 
DEMAND (22)  =  DTEMP (27) 
DEMAND (23)  =  DTEMP (14) 
DEMAND (24)  =  DTEMP (1) 
DEMAND(25)  =  DTEMP(24) 
DEMAND (26)  =  DTEMP (11) 
DEMAND (27)  =  DTEMP (21) 
DO  14  I  =  1,27,2 
D(I)  =  3 
DP (I)  =  1 
DPP(I)  =  2 
CONTINUE 
DO  15  1  =  2,26,2 
D(I)  =  1 
DP (I)  =  3 
DPP (I)  =  2 
CONTINUE 


ELSEIF ( (II ,EQ.KK) .AND. (II .GT. JJ) )THEN 


DEMAND (1) 
DEMAND ( 2 ) 
DEMAND ( 3 ) 
DEMAND (4) 
DEMAND  ( 5 ) 
DEMAND (6) 
DEMAND (7) 
DEMAND ( 8 ) 
DEMAND (9) 
DEMAND (10) 
DEMAND (11) 
DEMAND (12) 
DEMAND (13) 
DEMAND (14) 
DEMAND (15) 
DEMAND(16) 
DEMAND (17) 
DEMAND (18) 
DEMAND (19) 
DEMAND(20) 


=  DTEMP (21) 

=  DTEMP (12) 

=  DTEMP (20) 

=  DTEMP (24) 

=  DTEMP (11) 

=  DTEMP (3) 

=  DTEMP  (19) 

=  DTEMP (15) 

=  DTEMP (23) 

=  DTEMP (2) 

=  DTEMP (10) 
=  DTEMP (27) 
=  DTEMP (14) 
=  DTEMP (1) 

=  DTEMP (6) 

=  DTEMP (22) 
=  DTEMP (18) 
=  DTEMP(26) 
=  DTEMP(5) 

=  DTEMP  (13) 


DEMAND (21)  =  DTEMP(17) 

DEMAND (22)  =  DTEMP(4) 

DEMAND  (23)  =  DTEMPO) 

DEMAND (24)  =  DTEMP(25) 

DEMAND(25)  =  DTEMP(8) 

DEMAND(26)  =  DTEMP(16) 

DEMAND (27)  =  DTEMP(7) 

DO  16  I  =  1,13,2 
D(I)  =  2 
DP (I)  =  1 
DPP (I)  =  3 

16  CONTINUE 

DO  17  I  =  2,14,2 
D(I)  =  2 
DP (I)  =  3 
DPP (I)  =  1 

17  CONTINUE 

DO  116  I  =  15,27,2 
D(I)  =  2 
DP ( I )  =3 
DPP(I)  =  1 

116  CONTINUE 

DO  117  1  =  16,26,2 
D(I)  =  2 
DP (I)  =  1 
DPP (I)  =  3 

117  CONTINUE 

ELSEIF( (JJ.EQ.KK) .AND. (II .GT.JJ) )  THEN 

DEMAND(l)  =  DTEMPO) 

DEMAND(2)  =  DTEMP(21) 

DEMAND (3)  =  DTEMP(8) 

DEMAND(4)  =  DTEMP(20) 

DEMAND (5)  =  DTEMP(7) 

DEMAND (6)  =  DTEMP(19) 

DEMAND  (7)  =  DTEMP(6) 

DEMAND (8)  =  DTEMP(12) 

DEMANDO)  =  DTEMPdS) 

DEMAND(IO)  =  DTEMP(24) 

DEMAND(ll)  =  DTEMP(5) 

DEMAND(12)  =  DTEMP(ll) 

DEMAND (13)  =  DTEMP(17) 

DEMAND (14)  =  DTEMP(23) 

DEMAND (15)  =  DTEMP(4) 

DEMAND(16)  =  DTEMP(IO) 

DEMAND (17)  =  DTEMP(16) 

DEMAND(18)  =  DTEMP(22) 

DEMAND(19)  =  DTEMP(3) 

DEMAND(20)  =  DTEMP(15) 

DEMAND(21)  =  DTEMP(2) 

DEMAND (22)  =  DTEMP(27) 

DEMAND(23)  =  DTEMP(14) 

DEMAND(24)  =  DTEMP(l) 

DEMAND(25)  =  DTEMP(26) 

DEMAND (26)  =  DTEMP(13) 

DEMAND (27)  =  DTEMP(25) 

DO  18  I  =  1,27,2 
D(I)  =  3 
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DP (I)  =  2 
DPP(I)  =  1 

18  CONTINUE 

DO  19  I  =  2,26,2 
D(I)  =  2 
DP ( I )  =3 
DPP(I)  =  1 

19  CONTINUE 

ELSEIF{ (JJ.EQ.KK) .AND. (II.LT.JJ) )  THEN 

DEMAND(l)  =  DTEMP(25) 

DEMAND ( 2 )  =  DTEMP (16) 

DEMAND (3)  =  DTEMP (22) 

DEMAND (4)  =  DTEMP (13) 

DEMAND (5)  =  DTEMP (26) 

DEMAND (6)  =  DTEMP (7) 

DEMAND (7)  =  DTEMP (19) 

DEMAND (8)  =  DTEMP (17) 

DEMAND (9)  =  DTEMP (23) 

DEMAND (10)  =  DTEMP (4) 

DEMAND (11)  =  DTEMP (10) 

DEMAND (12)  =  DTEMP (27) 

DEMAND (13)  =  DTEMP (14) 

DEMAND (14)  =  DTEMP (1) 

DEMAND (15)  =  DTEMP (8) 

DEMAND(16)  =  DTEMP(20) 

DEMAND (17)  =  DTEMP (18) 

DEMAND(18)  =  DTEMP(24) 

DEMAND (19)  =  DTEMP (5) 

DEMAND (20)  =  DTEMP (11) 

DEMAND(21)  =  DTEMP(15) 

DEMAND (22)  =  DTEMP (2) 

DEMAND (23)  =  DTEMP (9) 

DEMAND (24)  =  DTEMP (21) 

DEMAND(25)  =  DTEMP(6) 

DEMAND(26)  =  DTEMP(12) 

DEMAND(27)  =  DTEMP(3) 

DO  20  I  =  1,13,2 
D(I)  =  1 
DP (I)  =  2 
DPP(I)  =  3 

20  CONTINUE 

DO  21  I  =  2,14,2 
D(I)  =  1 
DP (I)  =3 
DPP (I)  =  2 

21  CONTINUE 

DO  120  I  =  15,27,2 
D(I)  =  1 
DP (I)  =  3 
DPP (I)  =  2 

120  CONTINUE 

DO  121  I  =  16,26,2 
D(I)  =  1 
DP (I)  =2 
DPP(I)  =  3 

121  CONTINUE 
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ELSEIF( (II.LT.JJ) .AND. (JJ.LT.KK) )  THEN 

DEMAND(l)  =  DTEMP{2B) 

DEMAND (2)  =  DTEMP(16) 

DEMAND (3)  =  DTEMP(7) 

DEMAND (4)  =  DTEMP(22) 

DEMAND (5)  =  DTEMP(26) 

DEMAND (6)  =  DTEMP(13) 

DEMAND (7)  =  DTEMP(17) 

DEMAND (8)  =  DTEMP(4) 

DEMAND (9)  =  DTEMP(8) 

DEMAND(IO)  =  DTEMP(19) 

DEMAND(ll)  =  DTEMP(23) 

DEMAND(12)  =  DTEMP(IO) 

DEMAND (13)  =  DTEMP(27) 

DEMAND(14)  =  DTEMP(14) 

DEMAND (15)  =  DTEMP(l) 

DEMAND(16)  =  DTEMPdS) 

DEMAND (17)  =  DTEMP(5) 

DEMANDdS)  =  DTEMP(9) 

DEMAND(19)  =  DTEMP(20) 

DEMAND (20)  =  DTEMP(24) 

DEMAND(21)  =  DTEMPdl) 

DEMAND  (22)  =  DTEMPdS) 

DEMAND (23)  =  DTEMP(2) 

DEMAND (24)  =  DTEMP(6) 

DEMAND(25)  =  DTEMP(21) 

DEMAND (26)  =  DTEMP(12) 

DEMAND (27)  =  DTEMP(3) 

DO  22  I  =  1,27 
D(I)  =  1 
DP (I)  =  2 
DPP (I)  =  3 
22  CONTINUE 


ELSEIF( (II .LT 

DEMANDd) 
DEMAND (2) 
DEMAND (3) 
DEMAND (4) 
DEMAND (5) 
DEMAND (6) 
DEMAND (7) 
DEMAND(8) 
DEMAND (9) 
DEMAND (10) 
DEMAND (11) 
DEMAND (12) 
DEMAND (13) 
DEMAND (14) 
DEMANDdS) 
DEMANDdS) 
DEMAND (17) 
DEMANDdS) 
DEMAND (19) 
DEMAND(20) 
DEMAND (21) 
DEMAND (22) 
DEMAND(23) 


.KK) .AND. (KK 

=  DTEMP(25) 

=  DTEMP(22) 

=  DTEMP(19) 

=  DTEMP(16) 

=  DTEMP(26) 

=  DTEMPd3) 

=  DTEMP(23) 

=  DTEMPdO) 

=  DTEMP(20) 

=  DTEMP(7) 

=  DTEMP(17) 
=  DTEMP(4) 

=  DTEMP(27) 
=  DTEMP(14) 
=  DTEMPd) 

=  DTEMP(24) 
=  DTEMPdl) 
=  DTEMP(21) 
=  DTEMP(8) 

=  DTEMPdS) 
=  DTEMP(5) 

=  DTEMPdS) 
=  DTEMP(2) 


LT.JJ) )  THEN 
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DEMAND(24)  =  DTEKP(12) 

DE>1AND(25)  =  DTEMPO) 

DEMAND (26)  =  DTEMP(6) 

DE>1AND(27)  =  DTEMPO) 

DO  23  I  =  i,:;7 
D(I)  =  1 
DP ( I )  =  3 
DPP (I)  =  2 

23  CONTINUE 

ELSEIF( (JJ.UT.II) .AND. (II.LT.KK) )  THEN 

DEMAND(l)  =  DTEMP(21) 

DEMAND (2)  =  DTEMP{12) 

DEMAND (3)  DTEMP ( 3 ) 

DEMAND(4)  JTEMP(20) 

DEMAND (5)  -  DTEMP (24) 

DEMAND ( 6 )  =  DTEMP (11) 

DEMAND (7)  =  DTEMP (15) 

DEMAND ( 8 )  =  DTEMP ( 2 ) 

DEMAND (9)  =  DTEMP (6) 

DEMAND  (10)  DTEMP  (19) 

DEMAND(ll)  =  DTEMP (23) 

DEMAND (12)  =  DTEMP (10) 

DEMAND (13)  =  DTEMP (27) 

DEMAND (14)  =  DTEMP (14) 

DEMAND (15)  =  DTEMP(l) 

DEMAND (16)  =  DTEMP (18) 

DEMAND (17)  =  DTEMP (5) 

DEMAND (18)  =  DTEMP (9) 

DEMAND (19)  =  DTEMP (22) 

DEMAND(20)  =  DTEMP(26) 

DEMAND (21)  =  DTEMP (13) 

DEMAND (22)  =  DTEMP (17) 

DEMAND (23)  =  DTEMP (4) 

DEMAND (24)  =  DTEMP (8) 

DEMAND(25)  =  DTEMP(25) 

DEMAND(26)  =  DTEMP(16) 

DEMAND (27)  =  DTEMP (7) 

DO  24  I  =  1,27 
D(I)  =  2 
DP ( I )  =  1 
DPP(I)  =  3 

24  CONTINUE 

ELSEIF( (JJ.LT.KK) .AND. (KK.LT.II) )  THEN 


DEMAND(l) 
DEMAND (2) 
DEMAND (3) 
DEMAND (4) 
DEMAND(5) 
DEMAND (6) 
DEMAND (7) 
DEMAND(8) 
DEMAND (9) 
DEMAND (10) 
DEMAND (11) 
DEMAND (12) 
DEMAND (13) 


=  DTEMP (21) 

=  DTEMP (20) 

=  DTEMP (19) 

=  DTEMP (12) 

=  DTEMP (24) 

=  DTEMP (11) 

=  DTEMP (23) 

=  DTEMP(IO) 

=  DTEMP (22) 

=  DTEMPO) 

=  DTEMP (15) 
=  DTEMP (2) 

=  DTEMP (27) 
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DEKANDd-l)  =  DTEMP(14) 

DEa4Arro{i5)  =  dtemp(I) 

DEMAND(16)  =  DTEMP(26) 

DEMAND (17)  =  DTEMP(13) 

DEMAND (18)  -  DTEMP(25) 

DEMAND (19)  =  DTEMP(6) 

DEMAND(20)  =  DTEMPdS) 

DEMAND (21)  =  DTEMP(5) 

DEMAND (22)  =  DTEMP(17) 

DEMAND (23)  =  DTEMP(4) 

DEMAND(24)  =  DTEMP(16) 

DEMAND  (25)  =  DTEMPO) 

DEMAND (26)  =  DTEMP(8) 

DEMAND (27)  =  DTEMP(7) 

DO  25  I  =  1,27 
Dd)  =  2 
DPd)  =3 
DPP (I)  =  1 

25  CONTINUE 

ELSEIF( (KK.LT.II) .AND. (II.LT.JJ) )  THEN 

DEMANDd)  =  DTEMP(9) 

DEMAND (2)  =  DTEWr(6) 

DEMANDd)  =  DTEMi  (J) 

DEMAND (4)  =  DTEMP(8) 

DEMAND  (5)  -  DTEMPdd) 

DEMAND (6)  =  DTEMP(5) 

DEMAND(7)  =  DTEMP(15) 

DEMAND (8)  =  DTEMP(2) 

DEMANDd)  =  DTEMPd2) 

DEMAND (10)  =  DTEMP(7) 

DEMAND  (11)  =  DTEMPd?) 

DEMAND (12)  =  DTEMP(4) 

DEMAND (13)  =  DTEMP(27) 

DEMAND (14)  =  DTEMP(14) 

DEMAND  (15)  =  DTEMPd) 

DEMAND(16)  =  DTEMP(24) 

DEMAND  (17)  =  DTEMPdl) 

DEMANDdS)  =  DTEMP(21) 

DEMAND(19)  =  DTEMPd6) 

DEMAND(20)  =  DTEMP(26) 

DEMAND(21)  =  DTEMP(13) 

DEMAND (22)  =  DTEMP(23) 

DEMAND (23)  =  DTEMP(IO) 

DEMAND (24)  =  DTEMP(20) 

DEMAND(25)  ^  DTEMP(25) 

DEMAND (26)  =  DTEMP(22) 

DEMAND (27)  =  DTEMP(19) 

DO  26  I  =  1,27 
D(I)  =  3 
DP ( I )  =  1 
DPP (I)  =  2 

26  CONTINUE 

ELSEIF( (KK.LT.JJ) .AND. (JJ.LT.II) )  THEN 

DEMANDd)  =  DTEMPd) 

DEMAND(2)  =  DTEMP(8) 

DEMAND(3)  =  DTEMP(7) 
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DEMAND!4 ) 

=  DTEMP!6) 

DEMAND ! 5 ) 

-  DTEMP( 

18) 

DEMAND !6) 

=  DTEMP(5) 

DEMAND  17) 

=  DTEMP(17) 

DEMAND18) 

=  DTEMP(4) 

DEMAND !9) 

=  DTEMP( 

3) 

DEMAND  110) 

=  DTEMP 

(15) 

DEMAND  111) 

=  DTEMP 

(2) 

DEMAND! 12) 

=  DTEMP 

(27) 

DEMAND  113) 

=  DTEMP 

(14) 

DEMAND! 14) 

=  DTEMP 

(1) 

DEMAND! 15) 

=  DTEMP 

(26) 

DEMAND ! 1 6 ) 

=  DTEMP 

(13) 

DEMAND! 17) 

=  DTEMP 

(16) 

DEMAND ! 1 8 ) 

=  DTEMP 

(25) 

DEMAND! 19) 

=  DTEMP 

(12) 

DEMAND !20) 

=  DTEMP 

(24) 

DEMAND ! 2 1 ) 

=  DTEMP 

(11) 

DEMAND! 22) 

=  DTEMP 

(23) 

DEMAND! 23) 

=  DTEMP 

(10) 

DEMAND (24) 

=  DTEMP 

(22) 

DEMAND (25) 

-  DTEMP 

(21) 

DEMAND (26) 

=  DTEMP 

(20) 

DEMAND (27) 

=  DTEMP 

(19) 

DO  27  !  = 

1,27 

D!!)  = 

3 

DP  !  ! )  = 

2 

DPP(!) 

=  1 

CONTINUE 

ENDIF 

BEGIN  TO  FILL 

DEMAND  NODES 

DO  45  I  =  1,27 

IF(DEMAND(I) .E0.0)GOTO45 
IF(S(D(I))  .GE.  DEMAND(I))  THEN 
X(D(I) , I)  =  DEMAND! I) 

S(D(I) )  =  S(D(I) )  -  X(D(I) ,1) 
DEMAND ( I )  =0 
GOTO  45 

ELSEIF  (S(D(I))  .GT.  0)  THEN 
X(D(I) , I) =  S(D(I) ) 

S(D(I) )  =  S(D(I) )-  X(D(I) ,1) 
DEMAND! I)  =  DEMAND!!)  -  X!D!I),I) 
IF!DEMAND!I)  . EQ .  0)  G0T045 
ENDIF 


IF  1S!DP!I))  .GE.  DEMAND!!))  THEN 
X!DP!!) , !)  =  DEMAND!!) 

S!DP!!))  =  S!DP!!))  -  X(DP(I).!) 
DEMAND ! ! )  =0 
GOTO  45 


ELSEIF  !S!DP!!))  .GT.  0)  THEN 
X!DP!!) ,!)=  S!DP!!) ) 

S!DP!!))  =  S!DP!!))  -  X!DP!!) 
DEMAND!!)  =  DEMAND!!)  -  XIDP! 
!F!DEMAND!!) .EQ.O)GOTO  45 
ENDIF 


,!) 

I)  ,  I) 


IF(S(DPP(I))  -GE.  DEMAND(I))  THEN 
XCDPPd)  .1)  =  DEMAND(I) 

S(DPP(I))  =  S(DPP(I))  -  X(DPP(I),I) 

DEMAND! I)  =  0 
GOTO  45 

ELSEIF  (S(DPP(1))  .GT.O)  THEN 
X(DPP(I) , 1) =  S{DPP(I) ) 

S(DPP(I))  =  S(DPP(I))  -  X(DPP(I),I) 
DEMAND(I)  =  DEMAND! I)  -  X!DPP!I),I) 

END  IF 

IF  !DEMAND!I) .NE.O)  PRINT* ERROR DEMAND ! I ), I 
45  CONTINUE 


REORDER  FLOWS  TO  BE  INTERPOLATED  TO  THE  FINE  GRID  PROBLEM 

IF!  !II.EQ.JJ)  .AND.  !JJ.EQ.KK) )  THEN 
DO  30  L  =  1,DIM 

X2H!L, 1, 1, 1)  =  X!L.25) 

X2H!L,2,1,1)  =  X!L.22) 

X2H!L,3,1,1)  =  X!L,16) 

X2H(L, 1,2, 1)  =  X!L,23) 

X2H !L, 2 ,2,1)  =  X !L, 13 ) 

X2H!L,3,2,1)  =  X!L,4) 

X2H(L, 1,3,1)  =  X (L. 17) 

X2H!L,2,3,1)  =  X(L,5) 

X2H!L,3,3, 1)  =  X(L,1) 

X2H!L, 1, 1,2)  =  X!L,24) 

X2H!L,2, 1,2)  =  X(L, 14) 

X2H!L,3,1,2)  =  X!L,6) 

X2H(L, 1,2,2)  =  X(L,15) 

X2H!L,2,2,2)  =  X(L,26) 

X2H!L,3,2,2)  =  X!L.19) 

X2H!L, 1,3,2)  =  X(L,7) 

X2H!L,2,3,2)  =  X(L,20) 

X2H!L,3,3,2)  =  X!L,10) 

X2H!L, 1,1,3)  =  X!L,18) 

X2H!L,2,1,3)  =  X!L,8) 

X2H!L,3,1,3)  =  X!L,2) 

X2H!L,1,2,3)  =  X!L,9) 

X2H!L,2,2,3)  =  X!L,21) 

X2K!L,3,2,3)  =  X(L,11) 

X2H{L,1,3,3)  =  X(L,3) 

X2H(L,2,3,3)  =  X(L,12) 

X2H!L,3,3,3)  =  X!L,27) 

30  CONTINUE 

ELSEIF!  !II .EQ.JJ)  .AND. !KK.GT.JJ) )  THEN 
DO  31  L  =  1,DIM 

X2H!L, 1,1,1)  =  X(L,24) 

X2H(L,2, 1, 1)  =  X(L, 15) 

X2H!L,3,1,1)  =  X(L,5) 

X2HIL, 1,2,1)  =  X!L,16) 

X2H(L,2,2,1)  =  X!L.26) 

X2H!L,3,2,1)  =  X(L,17) 

X2H!L, 1,3,1)  =  X{L,6) 

X2H!L,2,3,1)  =  X(L,18) 

X2H!L,3,3,1)  =  X(L,27) 

X2H!L, 1,1,2)  =  X!L,21) 

X2H!L,2,1,2)  =  X!L,11) 
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X2H(L,3, 1,2)  =  X{L,3) 

X2H(L, 1,2,2)  =  X(L, 12) 
X2H(L,2,2,2)  =  X(L,23) 
X2H(L,3,2,2)  =  X(L,13) 

X2H(L, 1,3,2)  =  X(L,4) 

X2H (L, 2, 3 ,2)  =  X (L, 14 ) 

X2H(L,3,3,2)  -  X(L,251 

X2H(L, 1,1,3)  =  X(L, 19) 

X2H(L,2, 1,3)  =  X(L,7) 
X2H(L,3,1,3)  =  X(L,1) 
X2H(L,1,2,3)  =  X(L,8) 
X2H(L,2,2,3)  =  X(L,20) 
X2H(L,3,2,3)  =  X(L,9) 
X2H(L,1,3,3)  =  X(L,2) 
X2H(L,2,3,3)  =  X(L,10) 
X2H(L,3,3,3)  =  X(L,22) 

31  CONTINUE 

ELSEIF( (II ,EQ.JJ) .AND. (JJ.GT.KK) )  THEN 
DO  32  L  =  1,DIM 

X2H (L, 1,1,1)  =  X(L, 14) 

X2H (L, 2 ,1,1)  =  X (L, 10) 

X2H(L,3, 1,1)  =  X(L,6) 

X2H(I.,  1,2,1)  =  X(L,  11) 

X2H(L,,2,2, 1)  =  X(L,5) 

X2H (L, 3 ,2,1)  =  X (L, 2) 

X2H(L, 1,3,1)  =  X(L.7) 
X2H(L,2,3,1)  =  X(L,3) 
X2H(L,3,3,1)  =  X(L,1) 

X2H(L,1,] ,2)  =  X(L,22) 
X2H(L,2,1,2)  =  X(L,19) 

X2H(L,3, 1,2)  =  X(L,15) 

X2H(L, 1,2,2)  =  X(L.20) 
X2H(L,2,2,2)  =  X(L,13) 
X2H(L,3,2,2)  =  X(L,8) 

X2H(L, 1,3,2)  =  X(L,16) 
X2H(L,2,3,2)  =  X(L,9) 
X2H(L,3,3,2)  =  X(L,4) 

X2H(L, 1,1,3)  =  X(L,27) 
X2H(L,2,1,3)  =  X(L,25) 

X2H(L,3, 1,3)  =  X(L,23) 

X2H(L, 1,2,3)  =  X(L,26) 
X2H(L,2,2,3)  =  X(L.21) 
X2H(L,3,2,3)  =  X(L,17) 

X2H(L, 1,3,3)  =  X(L.24) 
X2H(L,2,3,3)  =  X(L.18) 
X2H(L,3,3,3)  =  X(L.12) 

32  CONTINUE 

ELSEIF( (II .EQ.KK) .AND. (JJ.GT.II) )  THEN 
DO  33  L  =  1,DIM 

X2H(L, 1,1,1)  =  X(L,24) 
X2H(L,2,1,1)  =  X(L,15) 
X2H(L,3,1,1)  =  X(L,5) 

X2H(L, 1,2,1)  =  X(L,21) 
X2H(L,2,2,1)  =  X(L,11) 
X2H(L,3,2,1)  =  X(L,3) 

X2H(L, 1,3,1)  =  X(L,19) 
X2H(L,2,3,1)  =  X(L,7) 

X2H(L,3,3, 1)  =  X(L,1) 

X2H(L, 1, 1,2)  =  X(L,16) 
X2H(L,2,1,2)  =  X(L,26) 
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X2H(L,3,1,2)  =  X(L,17) 

X2H(L, 1,2,2)  =  X(L,12) 

X2H(L,2,2,2)  =  X(L,23) 

X2H(L,3,2,2)  =  X(L,13) 

X2H(L, 1,3,2)  =  X(L,8) 

X2H(L,2,3,2)  =  X(L,20) 
X2H(L,3,3,2)  =  X(L.9) 

X2H(L, 1,1,3)  =  X(L,6) 

X2H(L,2, 1,3)  =  X(L, 18) 

X2H(L,3, 1,3)  =  X(L.27) 
X2H(L.1,2.3)  =  X(L,4) 

X2H(L,2.2,3)  =  X(L,14) 

X2H(L,3,2,3)  =  X(L,25) 

X2H(L, 1,3,3)  =  X(L.2) 

X2H(L,2.3,3)  =  X(L.IO) 
X2H(L,3,3,3)  =  X(L,22) 

33  CONTINUE 

ELSEIFt (lI.EO.KK) .AND. (II.GT.JJ) ) THEN 
DO  34  L  =  1,DIM 

X2H(L, 1, 1, 1)  =  X(L. 14) 
X2H(L,2,1,1)  =  X(L.IO) 
X2H(L,3,1,1)  =  X(L.6) 

X2H(L, 1,2,1)  =  X(L.22) 
X2H(L,2,2,1)  =  X(L,19) 
X2H(L,3,2,1)  =  X(L,15) 

X2H(L, 1,3,1)  =  X(L,27) 
X2H(L,2,3,1)  =  X(L,25) 
X2H(L,3,3,1)  =  X(L,23) 
X2H(L,1,1,2)  =  X(L,11) 
X2H(L,2,1,2)  =  X<L,5) 

X2H(L,3,1,2)  =  X(L,2) 

X2H(L, 1,2,2)  =  X(L,20) 
X2H(L,2,2,2)  =  X(L,13) 
X2H(L,3,2,2)  =  X(L,8) 

X2H(L, 1,3,2)  =  X(L,26) 
X2H(L,2,3,2)  =  X(L,21) 
X2H(L,3,3,2)  =  X(L,17) 

X2H(L, 1,1,3)  =  X(L,7) 
X2H(L,2,1,3)  =  X(L,3) 
X2H(L,3,1,3)  =  X(L,1) 

X2H(L, 1,2,3)  =  X(L, 16) 
X2H(L,2,2,3)  =  X(L,9) 
X2H(L,3,2,3)  =  X(L.4) 
X2H(L,1,3,3;  =  X(L,24) 
X2H(L,2,3,3)  =  X(L.18) 
X2H(L,3,3,3)  =  X(L,12) 

34  CONTINUE 

ELSEIF( (JJ.EQ.KK) .AND. (II.GT.JJ) )  THEN 
DO  35  L  =  1,DIM 

X2H(L, 1,1,1)  =  X(L,24) 
X2H(L,2,1,1)  =  X(L,21) 
X2H{L,3,1,1)  =  X(L,19) 

X2H(L, 1,2, 1)  =  X(L, 15) 
X2H(L,2,2,1)  =  X(L.ll) 
X2H(L,3,2,1)  =  X(L.7) 

X2H(L, 1,3,1)  =  X(L,5) 
X2H(L,2,3,1)  =  X(L.3) 
X2H(L,3,3,1)  =  X{L,1) 
X2H(L,1,1,2)  =  X(L.16) 
X2H(L,2,1,2)  =  X(L,12) 
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X2H{L,3, 1,2) 
X2H(L, 1,2,2) 
X2H(L,2,2,2) 
X2H(L,3,2,2) 
X2H(L, 1,3,2) 
X2H(L,2,3,2) 
X2H(L,3,3,2) 
X2H(L,1.1,3) 
X2H(L,2, 1,3) 
X2H(L,3, 1,3) 
X2H(L, 1,2,3) 
X2H(L,2,2,3) 
X2H(L,3,2,3) 
X2H(L, 1,3,3) 
X2H(L,2,3,3) 
X2H(L,3,3,3) 
CONTINUE 


X(L,8) 

X(L.26) 

X(L.23) 

X(L,20) 

X  (L,  17) 

X(L, 13) 

X(L,9) 

X(L.6) 

X(L,4) 

X(L,2) 

X (L, 18) 
X(L. 14) 
X(L, 10) 
X(L,27) 
X(L,25) 
X(L,22) 


ELSEIF( (JJ.EQ.KK) .AND. (II.LT.JJ) ) 
DO  36  L  =  1,DIM 


X2H(L, 1, 1, 1) 

= 

X(L, 14) 

X2H(L,2,1,1) 

= 

X (L,22) 

X2H(L,3, 1, 1) 

= 

X(L,27) 

X2H(L, 1,2,1) 

= 

X (L, 10) 

X2H(L,2,2, 1) 

= 

X(L,  19) 

X2H(L,3,2,1) 

= 

X(L,25) 

X2H(L, 1,3,1) 

= 

X(L,6) 

X2H(L,2,3, 1) 

:= 

X(L,15) 

X2H(L,3,3,1) 

= 

X(L,23) 

X2H(L,1, 1,2) 

= 

X(L,11) 

X2H(L,2, 1,2) 

s 

X(L,20) 

X2H(L,3, 1,2) 

= 

X(L,26) 

X2H(L, 1,2,2) 

= 

X(L,4) 

X2H{L,2,2,2) 

= 

X(L,13) 

X2H(L,3,2,2) 

X(L,21) 

X2H(L, 1,3,2) 

X(L,2) 

X2H(L,2,3,2) 

X(L,8) 

X2H(L,3,3,2) 

= 

X(L,17) 

X2H(L, 1,1,3) 

s 

X(L,7) 

X2H(L,2, 1,3) 

=s 

X(L, 16) 

X2H(L,3, 1,3) 

= 

X(L,24) 

X2H(L,1,2,3) 

= 

X(L,3) 

X2H(L,2,2,3) 

= 

X(L,9) 

X2H(L,3,2,3) 

= 

X(L, 18) 

X2H(L,1,3,3) 

= 

X(L,1) 

X2H(L,2,3,3) 

= 

X(L,5) 

X2H(L,3,3,3) 

= 

X(L, 12) 

CONTINUE 


ELSEIF( (II.LT.JJ) .AND. (JJ.LT.KK) ) 


DO  37  L  =  1,DIM 
X2H(L, 1,1,1) 
X2H(L,2,1,1) 
X2H{L,3, 1, 1) 
X2H(L, 1,2,1) 
X2H(L,2,2,1) 
X2H(L,3,2, 1) 
X2H(L, 1,3,1) 
X2H(L,2,3, 1) 
X2H(L,3,3,1) 
X2H(L,1,1,2) 
X2H(L,2,1,2) 


X(L,15) 

X(L,23) 

X(L,27) 

X(L,8) 

X(L,17) 

X{L,24) 

X(L,3) 

X(L,9) 

X(L, 18) 

X(L, 12) 

X(L,21) 


THEN 


THEN 


X2H(L,3,1,2)  =  X(L,26) 

X2H(L, 1,2,2)  =  X(L,6) 
X2H(L,2,2,2)  =  X(L,14) 
X2H(L,3,2,2)  =  X(L,22) 

X2H(L, 1,3,2)  =  X(L,2) 
X2H(L,2,3,2)  =  X(L,7) 
X2H(L,3,3,2)  =  X(L,16) 

X2H(L, 1,1,3)  =  X{L,10) 
X2H(L,2,1,3)  =  X(L,19) 

X2H(L,3, 1,3)  =  X(L,25) 

X2H(L, 1,2,3)  =  X{L,4) 
X2H(L,2,2,3)  =  X(L,11) 
X2H(L,3,2,3)  =  X(L.20) 
X2H(L,1,3,3)  =  X(L,1) 
X2H(L,2,3,3)  =  X(L,5) 
X2H(L,3,3,3)  =  X(L,13) 

37  CONTINUE 

ELSEIF( (II.LT.KK) .AND. (KK.LT.JJ) )  THEN 
DO  38  L  =  1,DIM 

X2H(L, 1,1,1)  =  X(L,15) 

X2H(L, 2, 1, 1)  =  X(L,23) 

X2H(L,3, 1,1)  =  X(L,27) 

X2H(L, 1,2, 1)  =  X(L, 12) 

X2H(L,2,2, 1)  =  X(L,21) 

X2H(L,3,2, 1)  =  X(L,26) 

X2H(L, 1,3,1)  =  X(L,10) 
X2H(L,2,3,1)  =  X(L,19) 
X2n(L,3,3,l)  =  X(L,25) 

X2H(L, 1,1,2)  =  X(L,8) 
X2H(L,2,1,2)  =  X<L,17) 
X2H(L,3,1,2)  =  X(L,24) 

X2H(L, 1,2,2)  =  X(L,6) 

X2H(L,2,2,2)  =  X(L,14) 
X2H(L,3,2,2)  =  X(L,22) 

X2H(L, 1,3,2)  =  X(L,4) 

X2H(L,2,3,2)  =  X(L.ll) 
X2H(L,3,3,2)  =  X(L,20) 

X2H(L, 1,1,3)  =  X(L.3) 

X2H(L,2,1,3)  =  X(L,9) 

X2H(L,3,1,3)  =  X(L.18) 

X2H(L, 1,2,3)  =  X  .L,2) 

X2H(L,2,2,3)  =  X(L,7) 

X2H(L,3,2,3)  =  X(L,16) 

X2H(L, 1,3,3)  =  X(L,1) 

X2H(L,2,3.3)  =  X(L,5) 

X2H(L,3,3,3)  =  X(L,13) 

38  CONTINUE 

ELSEIF( (JJ.LT. II) .AND. (II.LT.KK) )  THEN 
DO  39  L  =  1,DIM 

X2H(L, 1, 1, 1)  =  X(L.15) 
X2H(L,2,1,1)  =  X(L,8) 

X2H(L,3,1,1)  =  X(L,3) 

X2H(L, 1,2,1)  =  X(L,23) 
X2H(L,2,2,1)  =  X(L.17) 
X2H(L,3,2,1)  =  X(L,9) 

X2H(L, 1,3,1)  =  X(L,27) 
X2H(L,2,3,1)  =  X(L,24) 
X2H(L,3,3,1)  =  X(L,18) 

X2H(L, 1,1,2)  =  X(L,12) 
X2H(L,2,1,2)  =  X(L,6) 
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X2H(L,3, 1,2)  =  X(L,2) 

X2H(L, 1,2,2)  =  X(L,21) 
X2H(L,2,2,2)  =  X(L,14) 
X2H(L,3,2,2)  =  XIL,"?) 

X2H(L, 1,3,2)  =  X(L,26) 
X2H(L,2,3,2)  =  X(L,22) 
X2H(L,3,3,2)  =  X(L, 16) 

X2H (L, 1,1,3)  =  X (L, 10) 

X2H(L,2, 1,3)  =  X(L,4) 

X2H(L,3,1,3)  =  X(L,1) 

X2H(L, 1,2,3)  =  X(L, 19) 
X2H(L,2,2,3)  =  X(L, 11) 
X2H(L,3,2,3)  =  X(L,5) 
X2H(L,1,3,3)  =  X(L,25) 
X2H(L,2,3,3)  =  X(L,20) 
X2H(L,3.3,3)  =  X(L,13) 

39  CONTINUE 

ELSEIF( (JJ.LT.KK) .AND. (KK.LT.II) )  THEN 
DO  40  L  =  1,DIM 

X2H(L, 1,1,1)  =  X(L, 15) 
X2H(L,2,1,1)  =  X(L,12) 

X2H{L, 3,1,1)  =  X(L, 10) 

X2H(L, 1,2,1)  =  X(L,23) 
X2H(L,2,2,1)  =  X(L,21) 

X2H(L,3,2, 1)  =  X(L.19) 

X2H(L, 1,3,1)  =  X(L,27) 

X2H(L,2,3, 1)  =  X(L,26) 

X2H(L, 3,3,1)  =  X(L.25) 

X2H(L, 1.1,2)  =  X(L,8) 
X2H(L,2,1,2)  =  X(L,6) 

X2H(L,3,1,2)  =  X(L,4) 
X2H(L,1,2,2)  =  X(L,17) 
X2H(L,2.2,2)  =  X(L,14) 
X2H(L,3,2,2)  =  X(L,11) 
X2H(L,1,3,2)  =  X(L,24) 
X2H(L,2,3,2)  =  X(L,22) 
X2H(L,3,3,2)  =  X(L,20) 

X2H(L, 1,1,3)  =  X(L,3) 

X2H(L,2,1,3)  =  X(L.2) 

X2H(L,3,1,3)  =  X(L,1) 

X2H(L, 1,2,3)  =  X(L.9) 

X2H(L,2,2,3)  =  X(L,7) 

X2H(L,3,2,3)  =  X(L,5) 

X2H(L,1,3,3)  =  X(L,18) 
X2H(L,2,3,3)  =  X(L,16) 
X2H(L,3,3,3)  =  X(L,13) 

40  CONTINUE 

ELSEIF( (KK.LT.II) .AND. (II.LT.JJ) )  THEN 
DO  41  L  =  1,DIM 

X2H(L, 1,1,1)  =  X(L,15) 

X2H(L,2, 1, i)  =  X(L.8) 

X2H(L,3,1,1)  =  X(L,3) 

X2H(L, 1,2,1)  =  X(L,12) 
X2H(L,2,2,1)  =  X(L,6) 

X2H(L,3,2,1)  =  X(L,2) 

X2H(L, 1,3,1)  =  X(L,10) 
X2H(L,2,3,1)  =  X(L,4) 

X2H(L,3,3,1)  =  X(L,1) 

X2H{L, 1,1,2)  =  X(L,23) 
X2H(L,2,1,2)  =  X(L,17) 
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X2H(L,3,1,2)  =  X(L,9) 

X2H{L, 1,2,2)  =  X(L,21) 

X2H(L,2,2,2)  =  X(L,14) 

X2H(L,3,2,2)  =  X(L.7) 

X2H(L,1,3,2)  =  X(L,19) 

X2H(L,2,3,2)  =  X(L,11) 

X2H(L,3,3,2)  =  X(L,5) 

X2H(L,1,1,3)  =  X{L,27) 

X2H(L,2, 1,3)  =  X(L.24) 

X2H(L,3,1,3)  =  X(L.18) 

X2H(L, 1,2,3)  =  X(L,26) 

X2H(L,2,2,3)  =  X{L,22) 

X2H(L,3,2.3)  =  X(L,16) 

X2H(L. 1,3,3)  =  X{L,25) 

X2H(L,2,3,3)  =  X(L.20) 

X2H(L,3,3,3)  =  X(L.13) 

41  CONTINUE 

ELSEIF( (KK.LT.JJ) .AND. (JJ.LT. II) )  THEN 
DO  42  L  =  1,DIM 

X2H(L, 1,1,1)  =  X(L,14) 

X2H(L, 2, 1,1)  =  X(L, 11) 

X2H (L, 3 ,1,1)  =  X (L, 9) 

X2H(L, 1,2,1)  =  X(L,8) 

X2H(L,2,2, 1)  =  X(L.6) 

X2H(L,3.2, 1)  =  X(L,4) 

X2H(L, 1,3, 1)  =  X(L,3) 

X2H(L,2,3, 1)  =  X(L.2) 

X2H(L,3,3, 1)  =  X(L,1) 

X2H(L, 1, 1,2)  =  X(L,23) 

X2H(L,2,1,2)  =  X(L,21) 

X2H(L,3,1,2)  =  X(L,19) 

X2H(L, 1,2,2)  =  X(L,16) 

X2H(L,2,2,2)  =  X(L,13) 

X2H(L,3,2,2)  =  X(L,10) 

X2H(L, 1,3,2)  =  X(L,17) 

X2H(L,2,3,2)  =  X(L,7) 

X2H<L,3,3,2)  =  X(L,5) 

X2H(L, 1, 1,3)  =  X(L,27) 

X2H(L,2,1,3)  =  X(L,26) 

X2H(L,3,1,3)  =  X(L.25) 

X2H(L, 1,2,3)  =  X(L,24) 

X2H(L,2,2,3)  =  X(L.22) 

X2H(L,3,2,3)  =  X(L,20) 

X2H(L, 1,3,3)  =  X(L,18) 

X2H(L,2,3,3)  =  X(L,15) 

X2H{L,3,3,3)  =  X(L.12) 

42  CONTINUE 
ENDIF 
RETURN 
END 

C  THE  SUBROUTINE  INTER?  INTERPOLATES  THE  FLOWS  FROM  THE  COARSE  GRID 
C  BACK  TO  THE  FINE  GRID. 

C  THE  ARRAY  DH  CONTAINS  THE  DEMANDS  FOR  THE  FINE  GRID  DEMAND  NODES. 
C  THE  ARRAYS  DHTEMP,  XHTEMP,AND  STEM?  ARE  TEMPORARY  ARRAYS  USED 
C  TO  HOLD  THE  DEMAND,  FLOW  AND  SUPPLY  DURING  THE  INTERPOLATION 

C  SCHEME.  THE  ARRAY  X2H  CONTAINS  THE  FLOWS  FOR  THE  COARSE  GRID 

C  DEMAND  NODES.  THE  ARRAY  XH  CONTAINS  THE  INTERPOLATED  FLOW  FOR 

C  THE  FINE  GRID. 


80 


SUBROUTINE  INTERP  ( DH ,  DHTEMP ,  XHTEMP ,  STE^■iP ,  X2H  ,  XH ,  RC  ,  R ,  DIM ) 
INTEGER  DIM, RC, R, L, LL 

REAL  XH (DIM, R, R, R) , DH (R, R, R) , DHTEMP ( DIM , DIM , DIM ) , XP2H (3 , 3,3,3), 
tX2H (DIM, RC, RC, RC) , STEMP(DIM) , XHTEMP ( DIM , 2 ‘RC , 2 *RC , 2 ‘RC ) 

+ , IM, JM,KM, IP, JP,KP, IMJM, IMKM, IMJP, IMKP , JMKM , JMKP , J PKM , J PKP , 

+ IPJM, IPJP, IPKM, IPKP, IMJMKM, IMJMKP, IMJPKM, IMJPKP, IPJMKM, IPJMKP, 
+IPJPKM, IPJPKP 


DO  1  K  =  1,RC 

DO  1  J  =  1,RC 
DO  1  I  =  1,RC 

DO  1  L  =  1,DIM 

XHTEMP(L,2*I,2*J,2*K)  X2H(L,I,J,K) 

1  CONTINUE 
DO  2  K=l,3 

DO  2  J  =  1,3 
DO  2  I  =  1,3 

DO  2  L  =  1 , 3 

2  XP2H(L, I,J,K)=0 

C  RECOMPUTE  DEMANDS  FOR  THE  DEMAND  NODES  ON  THE  FINE  GRID 

DO  5  K  =  1,RC 
DO  5  J  =  1,RC 
DO  5  I  =  1,RC 
DO  6  KK  =  1,3 
DO  6  JJ  =  1,3 
DO  6  II  =  1,3 

DHTEMP(II,JJ,KK)  =  0 

6  CONTINUE 

DO  7  LL  =  1,DIM 

STEMP(LL)  =  XHTEMP(LL,2*I,2*J,2*K) 


CONTINUE 

IM  = 

.5 

IP  = 

.5 

JM  = 

.5 

JP  = 

.5 

KM  = 

.5 

KP  = 

.  5 

IMJM 

= 

.25 

IMJP 

= 

.25 

IPJM 

= 

.25 

IPJP 

= 

.25 

IMKM 

= 

.25 

IMKP 

= 

.25 

IPKM 

= 

.25 

IPKP 

= 

.25 

JMKM 

= 

.25 

JMKP 

.25 

JPKM 

= 

.25 

JPKP 

= 

.25 

IMJMKM 

=  .125 

IMJPKM 

=  .125 

IMJMKP 

=  .125 

IMJPKP 

=  .125 

IPJMKM 

=  .125 

IPJMKP 

=  .125 

IPJPKM 

=  .125 
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IPJPKP 


.125 


THE  IF  STATEMENTS  CHECK  FOR  EDGE  AND  CORNER  DEMAND  NODES. 

IF  (I  .EQ.  1)  THEN 
IM  =  1 
IMJM  =  .5 
IMJP  =  . 5 
IMKM  =  .5 
IMKP  =  .5 
IMJPKM  =  .25 
IMJMKP  =  .25 
IMJPKP  =  .25 
IMJMKM  =  .25 

ELSEIFd  .EQ.RC)  THEN 
IP  =  1 
I P JM  =  . 5 
IPJP  =  .5 
IPKM  =  .5 
IPKP  =  .5 
IPJPKM  =  .25 
IPJMKP  =  .25 
IPJPKP  =  .25 
IPJMKM  =  .25 
ENDIF 

IF  (J  .EO.  1)  THEN 
JM  =  1 
IMJM  =  .5 
IPJM  =  .5 
JMKM  =  .5 
JMKP  =  .5 
IMJMKM  =  .25 
IPJMKP  =  .25 
IMJMKP  =  .25 
IPJMKM  =  .25 
ELSEIF  (J  .EQ.RC)  THEN 
JP  =  1 
IMJP  =  .5 
IPJP  =  .5 
JPKM  =  .5 
JPKP  =  .5 
IMJPKM  =  .25 
IMJPKP  =  .25 
IPJPKP  =  .25 
IPJPKM  =  .25 
ENDIF 

IF  (K  .EQ.  1)  THEN 
KM  =  1 
IMKM  =  .5 
IPKM  =  .5 
JMKM  =  .5 
JPKM  =  .5 
IMJPKM  =  .25 
IMJMKM  =  .25 
IPJPKM  =  .25 
IPJMKM  =  .25 
ELSEIF  (K  .EQ.RC)  THEN 
KP  =  1 
IMKP  =  .5 


IPKP  = 

.5 

JMKP  = 

.5 

JPKP  = 

.  5 

IMJPKP 

.  25 

IMJMKP 

= 

.25 

IPJPKP 

r 

.25 

IPJMKP 

= 

.25 

END  IF 


DHTEMP (2 
DHTEMP ( 1 
DHTEMP ( 3 
DHTEMP (2 
DHTEMP (2 
DHTEMP (2 
DHTEMP (2 


.2.2)  = 
.2.2)  = 
.2.2)  = 
,1,2)  = 
,3,2)  = 
,2,1)  = 
.2,3)  = 


DH ( 2 • I , 2 
IM'DH (2* 
IP*DH (2* 
JM*DH(2* 
JP*DH(2* 
KM*DH (2* 
KP«DH(2« 


IF  ( (I  .EQ. 1)  .AND.  (J  .EQ. 
IMJM  =  1 
IMJMKP  =  . 5 
IMJMKM  =  . 5 

ELSEIF  ((I  .EQ.l)  .AND.  (J 
IMJP  =  1 
IMJPKP  =  .5 
IMJPKM  =  .5 

ENDIF 

IF  {{I  .EQ.RC)  .AND.  (J  .EO 
IPJM  =  1 
IPJMKP  =  .5 
I P JMKM  =  . 5 

ELSEIF  ((I  .EQ.RC)  .AND.  (J 
IPJP  =  1 
IPJPKP  =  .5 
IPJPKM  =  .5 

ENDIF 

IF  ((I  .EQ.l)  .AND.  (K  .EQ. 
IMKM  =  1 
IMJMKM  =  .5 
IMJPKM  =  .5 

ELSEIF  ((I  .EQ.l)  .AND.  (K  . 
IMKP  =  1 
IMJMKP  =  .5 
IMJPKP  =  .5 

ENDIF 

IF  ((I  .EQ.RC)  .AND.  (K  .EQ. 
IPKM  =  1 
IPJMKM  =  .5 
IPJPKM  =  .5 

ELSEIF  ((I  .EQ.RC)  .AND.  (K 
IPKP  =  1 
IPJMKP  =  .5 
IPJPKP  =  .5 

ENDIF 

IF  ( ( J  .EQ.  1)  .AND.  (K  .EQ. 
JMKM  =  1 
IMJMKM  =  .5 
IPJMKM  =  .5 

ELSEIF  ( (J  .EQ.  1)  .and.  (K 
JMKP  =  1 
IMJMKP  =  ,  5 


*J,2*K) 

I-1,2*J,2»K) 

I+1,2*J,2*K) 

I,2*J-1,2*K) 

I,2*J+1,2*K) 

1,2*J,2»K-1) 

I,2*J,2*K+1) 

1 ) )  THEN 


.EQ.RC))  THEN 


1 ) )  THEN 


.EQ.RC))  THEN 


1 ) )  THEN 


EQ.RC))  THEN 


1 ) )  THEN 


.EQ.RC))  THEN 


1) )  THEN 


.EQ.RC))  THEN 
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IPJMKP  =  .5 
END  IF 

IF  ( (J  .EQ.RC)  .AND.  (K  . EQ .  1))  THEN 

JPKM  =  1 
IMJPKM  =  .5 
IPJPKM  =  .5 

ELSEIF  ( (J  .EQ.RC)  .AND.  (K  .EQ.RC))  THEN 
JPKP  =  1 
IMJPKP  =  .5 
IPJPKP  =  .5 
ENDIF 

DHTEMP (1, 1 , 2)  =  IMJM*DH(2*I-1,2*J-1,2*K) 
DHTEMP{1,3 ,2)  =  IMJP'DH ( 2  *  I - 1 , 2  * J  + 1 , 2 *K ) 
DHTEMP(1,2, 1)  =  IMKM*DH(2*I-1,2»J,2«K-1) 
DHTEMP(1,2,3)  =  IMKP«DH (2*1-1 , 2*J, 2«K+1 ) 

DHTEMP (3 . 1 , 2)  =  I PJM*DH ( 2 *  I  + 1 , 2 * J - 1 , 2 'K ) 

DHTEMP (3 , 3 , 2)  =  IPJP*DH ( 2  *  I  + 1 , 2 * J  + 1 , 2 ‘K ) 
DHTEMP(3,2,1)  =  IPKM*DH ( 2 *  I  + 1 , 2 * J , 2 ‘K - 1 ) 
DHTEMP(3,2,3)  =  I PKP*DH ( 2  *  I + 1 , 2  * J , 2 *K  + 1 ) 
DHTEMP(2, 1, 1)  =  JMKM*DH(2*I.2*J-1,2*K-1) 
DHTEMP(2,1,3)  =  JMKP*DH ( 2 *  1 , 2 * J - 1 , 2 *K+ 1 ) 
DHTEMP(2,3, 1)  =  JPKM*DH ( 2  *  1 , 2 * J  + 1 , 2 ‘K -  1 ) 
DHTEMP(2,3,3)  =  JPKP*DH ( 2  *  1 , 2 * J  +  1 , 2 ‘K  +  1 ) 

IF  ( (I.EQ.l) .AND. (J.EQ.l) .AND. (K.EQ.l) )  THEN 
IMJMKM  =  1 

ELSEIF  ( (I .EQ. 1) .AND. (J.EQ.l) .AND. (K. EQ.RC) )  THEN 
IMJMKP  =  1 

ELSEIF  ( (I.EQ.l) .AND. (J. EQ.RC) .AND. (K.EQ.l) )  THEN 
IMJPKM  =  1 

ELSEIF  ( (I.EQ.l) .AND. (J. EQ.RC) .AND. (K. EQ.RC) )  THEN 
IMJPKP  =  1 

ELSEIF  ( (I. EQ.RC) .AND. (J.EQ.l) .AND. (K.EQ.l) )  THEN 
IPJMKM  =  1 

ELSEIF  ( (I. EQ.RC) .AND. (J. EQ.RC) .AND. (K.EQ.l) )  THEN 
IPJPKM  =  1 

ELSEIF  ( (I. EQ.RC) .AND. (J. EQ.RC) .AND. (K. EQ.RC) )  THEN 
IPJPKP  =  1 

ELSEIF  ( (I .EQ.RC) .AND. (J.EQ. 1) .AND. (K. EQ.RC) )  THEN 
IPJMKP  =  1 
ENDIF 

DHTEMPd,  1, 1)  =IMJMKM*DH(2*I-1,2*J-1,2*K-1) 
DHTEMP(1,1,3) =IMJMKP*DH(2*I-1,2*J-1,2*K+1) 
DHTEMP(1,3,1) =IMJPKM*DH(2*I-1,2*J+1,2*K-1) 
DHTEMP(1,3,3) =IMJPKP*DH(2*I-1,2*J+1,2*K+1) 

DHTEMP (3 , 1, 1) =IPJMKM*DH(2*I+1,2*J-1,2*K-1) 

DHTEMP (3, 1,3) =IPJMKP*DH (2* 1+1,2* J-1,2*K+1) 
DHTEMP(3,3,1)=IPJPKM*DH(2*I+1,2*J+1,2*K-1) 
DHTEMP(3,3,3) =IPJPKP*DH(2*I+1,2*J+1,2*K+1) 

IT  =  2*1-1 
JT  =  2*J-1 
KT  =  2*K-1 

SOLVE  THE  LOCAL  3X27  PROBLEM 

CALL  LSOLVE(XP2H, DHTEMP, STEMP, DIM, IT, JT,KT) 

DO  25  L  =  1,DIM 
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C  ACCUMULATE  THE  FLOW  FOR  EACH  FINE  GRID  DEMAND  NODE. 

XH (L, 2*1 -I , -1 , 2*K-1 ) =XP2H (L, 1,1,1) +XH ( L , 2  *  I  -  1 , 2  * J -  1 , 2  * K - 1 ) 

XH (L, 2*1 , 2*J -  1 , 2*K- 1 ) =XP2H (L, 2, 1 , 1 ) +XH ( L , 2  *  I , 2 • J -  1 . 2*K- 1 ) 

XH  (L,  2*  I  f  1 , 2*J  -1 , 2*K  -  1  )  =XP2H  (L,  3,l,l)tXH(L,2*Ii-l,2-J-l,2*K-l) 
XH(L,2«I-1,2*J,2»K-1) =XP2H (L,1,2,1)*XH(L,2*I-1,2*J,2*K-1) 
XH(L,2*I,2»J,2*K-1)=XP2H(L,2,2,1) tXH ( L , 2 • I , 2 • J , 2*K- 1 i 
XH (L, 2*1 , 2«U , 2*K-1 ) -XP2H (L, 3,2,1)  *XH ( L , 2  *  1 + 1 , 2  * J , 2*K  1 ) 

XH  (L,  2*  I  -1 , 2*J^-1  ,  2*K  -  1  )  =XP2H  (L,  1,3,1)  +XH  (L,  2*1  -1 , 2  *  J  ♦  1 , 2  *K  -  1  ) 
XH(L,2*I,2*J  +  1,2*K-1)  XP2H  (L,  2, 3 , 1 )  tXH  (L,  2*1 , 2*Jtl  ,  2*K  - 1  ) 

XH  (L,  2*1  tl  ,  2*J-fl  ,  2*K  -1 )  =  XP2H(L,  3,3,1)  +XH  (L,  2*  W  1,2*J  +  1,2*K-1) 
XH(L, 2*1-1, 2*J-1,2*K)  =  XP2H(L,1,1,2)+XH(L, 2*1-1, 2*J-1,2*K) 

XH (L, 2*1 , 2*J-1 , 2*K)  =  XP2H (L, 2, 1 , 2) +XH{L, 2*1 , 2*J-1 , 2*K) 

XH(L. 2*1+1, 2*J-1,2*K)  =  XP2H (L, 3 , 1, 2) +XH(L, 2*1+1 , 2*J- 1 , 2*K) 

XH(L, 2*1-1, 2*J,2*K)  =  XP2H (L, 1 , 2, 2) +XH (L, 2*1 -1 , 2*J , 2*K) 

XH(L,2*I,2*J,2*K)  =  XP2H(L,2,2,2)+XH(L,2*I,2*J,2*K) 

XH(L, 2*1+1, 2*J,2*K)  =  XP2H(L,3,2,2) +XH(L,2*I+1,2*J,2*K) 

XH(L, 2*1-1, 2*J+1,2*K)  =  XP2H (L, 1 , 3 , 2) +XH (L, 2*1-1 , 2*J+1 , 2*K) 

XH(L, 2*1, 2*1+1, 2*K)  =  XP2H(L,2,3,2)+XH(L.2*I,2*J+1,2*K) 

XH{L, 2*1+1, 2*J+1,2*K)  =  XP2H(L,3.3,2) +XH(L,2*I+1,2*J+1,2*K) 

XH(L, 2*1-1, 2*J-1,2*K+1)  =  XP2H(L, 1 , 1 ,3 ) +XH (L, 2*1-1 , 2*J-1 , 2*K+1 ) 
XH(L,2*I,2*J-1,2*K+1)  =  XP2H (L, 2, 1 , 3 ) +XH (L, 2*1 , 2*J-1 , 2*K+ 1 ) 

XH{L,2*I  +  1,2*J-1,2*K+1)  =  XP2H (L, 3.1,3) +XH ( L , 2  *  I  +  1 , 2  * J -  1 , 2*K  +  1 ) 
XH(L, 2*1-1, 2*J,2*K+1)  =  XP2H (L, 1, 2, 3 ) +XH (L, 2*1 -1 , 2*J , 2*K+ 1 ) 

XH(L,2*I,2*J,2*K+1)  =  XP2H{L,2,2,3) +XH(L,2*I,2*J,2*K+1) 

XH(L, 2*1+1, 2*J,2*K+1)  =  XP2H(L,3,2,3)+XH(L,2*I+1,2*J,2*K+1) 

XH(L, 2*1-1, 2*J+1,2*K+1)  =  XP2H(L,1.3,3>+XH(L,2*I-1,2*J+1,2*K+1) 

XH (L, 2*1 , 2*J+1 , 2*K+1 )  =  XP2H(L,2,3,3)+XH(L,2*I,2*J+1,2*K+1) 

XH(L,2*I+1,2*J+1,2*K+1)=XP2H(L,3,3,3)+XH(L,2*I+1,2*J+1,2*K+1) 

25  CONTINUE 

5  CONTINUE 
RETURN 
END 

C  THE  SUBROUTINE  RELAX  IMPROVES  THE  SOLUTION  OBTAINED  BY  THE 

C  THREE  PREVIOUS  SUBROUTINES.  THE  SUBROUTINE  CHECKS  THE  FLOW  TO 

C  EACH  DEMAND  FOR  THE  CHEAPEST  SUPPLIER.  IF  THE  FLOW  IS  NOT  FROM 

C  THE  CHEAPEST  SUPPLIER,  A  DUMMY  PROBLEM  IS  CREATED  THAT  IS  OF  THE 

C  FORM  AS  THE  ORIGINAL  PROBLEM.  THE  ONLY  SUPPLY  AND  DEMAND  IS 

C  FROM  THE  MISDIRECTED  FLOW.  THE  ARRAY  STEMP  CONTAINS  ALL  THE 

C  MISDIRECTED  FLOW.  THE  ARRAY  XH  CONTAINS  THE  FLOW  FOR  EACH  FINE 

C  GRID  DEMAND  NODE.  THE  ARRAY  DPTEMP  CONTAINS  THE  MISFILLED  DEMAND 

C  NODE'S  DEMAND. 

C  THE  RELAXATION  TAKES  PLACE  ON  THE  GLOBAL  GRID  BY  SOLVING  SEVERAL 

C  LOCAL  3X27  PROBLEMS. 

SUBROUTINE  RELAX ( STEMP , XH , RC , R , DIM) 

INTEGER  DIM,R,RC,M,L,0,P,Q,S,C 

REAL  STEMP(DIM) , DPTEMP ( 3 , 3 , 3 ) , XH (DIM, R , R , R } .COST, 
+DTEMP(3,63,63,63) , XHH (3 , 63 , 63 , 63 ) .SUM, 

+XTEMP ( 3 , 3 , 3 , 3 ) 

CC=0 

DO  10  K=1,R 

DO  10  J  =  1,R 
DO  10  I  =  1,R 

DO  10  L  =  1,DIM 

DTEMP(L, I, J,K)  =  0 
XHH(L,I,J,K)  =  0 

10  CONTINUE 

DO  15  K=l,3 
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n  n 


DO  15  J=:l,3 
DO  15  1=1,3 

DPTEMPd,  J,K)  =0 
DO  15  L=1,DIM 

XTEMP(L, I,J,K) =0 

15  CONTINUE 


STEMP ( 1 )  =0 
STEMP(2)  =  0 
STEMP (3)  =  0 
COST  =  0 
SUM  =  0 

OPEN (UNIT  =  12,  FILE  =  ' FLOCOST  DATA  D‘, STATUS  =  ‘OLD') 

CHECK  THE  FLOW  TO  EACH  DEMAND  NODE  TO  VERIFY  IF  THE  CHEAPEST 
SUPPLIER  WAS  USED.  IF  NOT  CORRECT  THE  FLOW. 

DO  20  K  =  1,R 
DO  20  J  =  1,R 
DO  20  I  =  1,R 

DO  20  L  =  1,DIM 
M  =  MIN(I,J,K) 

IF(L  .EO.  1)  THEN 
IF(I.GT.M)THEN 

DTEMP(L, I, J.K)  =  XH(L,I,J,K)+DTEMP(L,I,J,K) 
STEMP(l)  =  XH(L,I,J,K)+  STEMP ( 1 ) 

XH(L,I,J.K)  =  0 

COST  =  COST  +  DTEMP(L, I, J,K) *I 
END  IF 

COST  =  COST  +  I*XH(L, I, J.K) 

ELSEIF  (L  .EQ.  2)  THEN 
IF(J.GT.M)THEN 

DTEMP(L, I, J.K)  =  XH(L.I.J.K)+DTEMP{L.I,J,K) 
STEMP(2)  =  XH(L.I.J.K)  +  STEMP(2) 
XH(L.I.J.K)  =  0 

COST  =  COST  +  DTEMP(L. I. J.K) *J 
END  IF 

COST  =  COST  +  J*XH(L. I. J.K) 

ELSE 

IF (K.GT.M)THEN 

DTEMP(L. I.J.K)  =  XH(L. I.J.K) +DTEMP(L. I. J.K) 
STEMP (3)  =  XH(L. I.J.K) +  STEMP (3) 

XH(L. I.J.K)  =  0 

COST  =  COST  +  DTEMP(L. I.J.K) *K 
ENDIF 

COST  =  COST  +  K*XH(L. I.J.K) 

ENDIF 

SUM  =  SUM  +  XH(L. I.J.K)  *  DTEMP (L. I . J.K) 

20  CONTINUE 

WRITE(12. 13)  COST. SUM. STEMP (1) .STEMP(2) .STEMP(3) 

12  FORMAT (2X. 12. 2X. 12. 2X. 12 . 2X. 12 . 2X . F9 . 2 . 2X. F9 . 2 ) 

13  FORMAT ( 2X , F12 , 2 . 2X . F12 . 2 . 2X . F8 . 2 , 2X . F8 . 2 . 2X . F8 , 2 ) 

DO  30  K  =  l.RC 
DO  30  J  =  l.RC 
DO  30  I  =1.RC 

DO  30  L  =  l.DIM 

DPTEMPd.  1. 1)  =  DTEMP(L, 2*1-1. 2*J-1.2*K-1) 
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DTEMP (L, 2*1 -1 , 2*J-1 , 2»K-1 )  =  0 

DPTEMP(2,1,1)  =  DTEMP(L,2*I,2*J-1,2*K-1) 
DTEMP(L,2*I,2*J-1,2*K-1)  =  0 

DPTEMP(3,1,1)  =  DTEMP(L, 2*1+1, 2*J-1,2*K-1) 

DTEMP(L, 2*1+1, 2*J-1,2*K-1)  =  0 

DPTEMP (1 , 2 , 1 )  =  DTEMP(L, 2*1-1, 2*J,2*K-1) 

DTEMP(L, 2*1-1, 2*J,2*K-1)  =  0 

DPTEMP(2,2, 1)  =  DTEMP(L,2*I,2*J,2*K-1) 
DTEMP(L,2*I,2*J,2*K-1)  =  0 
DPTEMP(3,2,1)  =  DTEMP(L, 2*1+1, 2*J,2*K-1) 

DTEMP(L, 2*1+1, 2*J,2*K-1)  =  0 

DPTEMP(1,3,1)  =  DTEMP(L, 2*1-1, 2*J+1,2*K-1) 

DTEMP(L, 2*1-1, 2*J+1,2*K-1)  =  0 
DPTEMP(2,3,1)  =  DTEMP(L,2*I,2*J+1,2*K-1) 
DTEMP(L,2*I,2*J+1,2*K-1)  =  0 

DPTEMP(3,3,1)  =  DTEMP(L, 2*1+1, 2*J+1,2*K-1) 

DTEMP(L, 2*1+1, 2*J+1,2*K-1)  =  0 
DPTEMPd,  1,2)  =  DTEMP(L, 2*1-1, 2*J-1,2*K) 

DTEMP(L, 2*1-1, 2*J-1,2*K)  =  0 
DPTEWP(2, 1,2)  =  DTEMP(L,2*I,2*J-1,2*K) 
DTEMP(L,2*I,2*J-1,2*K)  =  0 
DPTEMP(3,  1,2)  =  DTEl>iP(L, 2*1  +  1, 2*J-1,2*K) 

DTEMP(L, 2*1+1, 2*J-1.2*K)  =  0 
DPTEMPd, 2, 2)  =  DTEMP(L, 2*1-1, 2*J,2*K) 

DTEMP(L, 2*1-1, 2*J,2*K)  =  0 
DPTEMP(2,2,2)  =  DTEMP (L, 2*1 , 2*J , 2 *K) 
DTEMP(L,2*I,2*J.2*K)  =  0 
DPTEMP(3 ,2,2)  =  DTEMP (L , 2 *  I +  1 , 2* J , 2 *K) 

DTEMP{L, 2*1+1, 2*J,2*K)  =  0 

DPTEMPd, 3, 2)  =  DTEMP(L, 2*1-1, 2*J+1,2*K) 

DTEMP(L, 2*1-1, 2*J+1,2*K)  =  0 

DPTEMP (2, 3, 2)  =  DTEMP (L. 2*1 , 2 • J+1 , 2*K) 

DTEMP(L,2*I,2*J+1,2*K)  =  0 

DPTEMPd, 3, 2)  =  DTEMP(L. 2*1  +  1, 2*J+1,2*K) 

DTEMP(L, 2*1+1, 2*J+1.2*K)  =  0 

DPTEMPd, 1,3)  =  DTEMP(L, 2*1-1, 2*J-1,2*K+1) 

DTEMP(L, 2*1-1, 2*J-1,2*K+1)  =  0 
DPTEMP(2, 1,3)  =  DTEMP(L,2*I,2*J-1,2*K+1) 
DTEMP(L,2*I,2*J-1,2*K+1)  =  0 
DPTEMP(3, 1,3)  =  DTEMP(L, 2*1+1, 2*J-1,2*K+1) 

DTEMP(L, 2*1+1, 2*J-1,2*K+1)  =  0 
DPTEMPd, 2, 3)  =  DTEMP(L,2*I-  ,2*J,2*K+1) 

DTEMP(L, 2*1-1, 2''J,2*K+1)  =  0 
DPTEMP(2,2,3)  =  Dl'EMP  (L.  2*1 , 2*J  ,  2*K+1 ) 
DTEMP(L,2*I,2*J,2*K+1)  =  0 
DPTEMP(3,2,3)  =  DTEMP (L, 2*1+1 , 2*J , 2*K+1 ) 

DTEMP (L, 2*1+1, 2*J,2*K+1)  =  0 

DPTEMPd, 3, 3)  =  DTEMP(L, 2*1-1, 2*J  +  1,2*K+1) 

DTEMP(L, 2*1-1, 2*J+1.2*K+1)  =  0 
DPTEMP(2,3,3)  =  DTEMP (L, 2*1 , 2*J+1 , 2*K+1 ) 
DTEMP(L,2*I,2*J+1,2*K+1)  =  0 
DPTEMP(3,3,3)  =  DTEMP (L, 2 *1+1 , 2*J+1 , 2*K+1 ) 

DTEMP(L, 2*1+1, 2*J+1.2*K+1)  =  0 
ITEMP  ^  2*1-1 
JTEMP  =  2*J-1 
KTEMP  =  2*K-1 

C  SOLVE  THE  LOCAL  3X27  DUMMY  PROBLEM  TO  CORRECT  THE  FLOW. 

CALL  LSOLVE(XTEMP, DPTEMP, STEMP, DIM, ITEMP, JTEMP, KTEMP) 
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C  ACCUMULATE  THE  FLOW  FOR  EACH  FINE  GRID  DEMAND  NODE. 

DO  140  C  =  1,3 

XHH ( C , 2  » I  - 1 , 2 • J -  1 , 2 ‘K - 1 ) =XTEMP (C , 1 , 1 , 1 ) +XHH (C,2«I-1,2»J-1,2*K-1) 
XHH(C,2*I,2*J-1,2*K-1)  =  XTEMP (C , 2 , 1 , 1 WXHH (C , 2» I , 2* J -  1 , 2*K- 1 ) 
XHH (C,2»I+1,2*J-1,2*K-1) =XTEMP (C, 3 , 1 , 1 ) +XHH (C,2*I+1,2»J-1,2*K-1) 
XHH (C, 2*1 -1 , 2*J , 2*K-1 )  =  XTEMP (C , 1 , 2 , 1 ) +XHH (C , 2* I  -  1 , 2  * J , 2 *K - 1 ) 
XHH(C,2*I,2*J,2*K-1)  =  XTEMP (C , 2 , 2 , 1 ) +XHH (C , 2  *  1 , 2  * J , 2 ‘K -  1 ) 

XHH (C, 2*1  +  1 , 2*J , 2*K-1 )  =  XTEMP (C , 3 , 2 . 1 ) +XHH ( C , 2  *  I  +  1 , 2  * J  ,  2 *K -  1 ) 
XHH (C,2*I-1,2*J+1,2*K-1) =XTEMP (C, 1 , 3 , 1 ) +XHH (C, 2*1-1, 2*Jtl,2*K-l) 
XHH(C,2*I,2*J+1,2*K-1)=XTEMP(C,2,3,1) +XHH(C,2*I,2*J+1,2*K-1) 

XHH (C, 2*1+1 , 2*J+1 , 2*K-1) =XTEMP(C, 3 ,3,1) +XHH (C,2*I+1,2*J+1,2*K-1) 
XHH (C, 2*1 -1 , 2*J-1 , 2*K)  =  XTEMP (C , 1 , 1 , 2 ) +XHH (C , 2 *  I - 1 , 2  * J - 1 , 2 *K ) 
XHH(C,2*I,2*J-1,2*K)  =  XTEMP (C, 2 , 1 , 2 ) +XHH (C, 2  *  1 , 2 *J - 1 , 2 *K ) 

XHH (C, 2*1 +1 , 2*J-1 , 2*K)  =  XTEMP (C , 3 , 1 , 2 ) +XHH (C , 2 *  I  + 1 , 2 * J - 1 , 2 *K ) 
XHH(C, 2*1-1, 2*J,2*K)  =  XTEMP (C , 1 , 2 , 2 ) +XHH (C , 2* I - 1 , 2  * J  ,  2 *K ) 

XHH (C, 2*1 , 2*J , 2*K)  =  XTEMP (C , 2 , 2 , 2 ) +XHH (C , 2 *  1 , 2*J , 2*K ) 

XHH (C, 2* I  +  l , 2*J , 2*K)  =  XTEMP (C , 3 , 2 , 2 ) +XHH (C , 2 *  I  + 1 , 2 * J , 2 *K ) 

XHH(C, 2*1-1, 2*J  +  1,2*K)  =  XTEMP (C, 1 , 3 , 2 ) +XHH (C , 2* I - 1 , 2  * J  + 1 , 2 *K ) 
XHH(C,2*I,2*J  +  1,2*K)  =  XTEMP (C , 2 , 3 , 2 ) +XHH (C , 2  *  I , 2  * J  + 1 , 2 *K) 

XHH(C, 2*1  +  1 , 2*J  +  1 , 2*K)  =  XTEMP (C , 3 , 3 , 2 ) +XHH (C , 2  *  I  +  1 , 2  * J  +  1 , 2 *K ) 
XHH (C, 2* I  -  1 , 2*J- 1 , 2*K  +  1) =XTEMP(C, 1,1,3) +XHH (C , 2 *  I -  1 , 2 * J -  1 , 2*K+1 ) 
XHH(C,2*I,2*J-1,2*K  +  1)  =  XTEMP (C , 2 , 1 , 3 ) +XHH (C , 2  *  1 , 2  * J -  1 , 2 *K  +  1 ) 
XHH (C,2*I+1.2*J-1,2*K+1) =XTEMP(C, 3 ,1,3) +XHH (C,2*I+1,2*J-1,2*K+1) 
XHH(C, 2*1-1, 2*J,2*K  +  1)  =  XTEMP (C, 1 , 2 , 3 ) +XHH (C , 2* I - 1 , 2  * J , 2 *K+ 1 ) 
XHH(C,2*I,2*J,2*K+1)  =  XTEMP (C, 2 , 2 , 3 ) +XHH (C , 2* I , 2* J , 2*K+ 1 ) 

XHH(C, 2*1  +  1, 2*J,2*K+1)  =  XTEMP (C, 3 , 2 , 3 ) +XHH (C , 2  * I  +  l , 2*J , 2 *K  +  1 ) 
XHH(C,2*I-l,2*J+l,2*K+l)=XTEMPvC, 1,3.3) +XHH(C,2*I-1,2*J+1,2*K+1) 
XHH(C,2*I,2*J  +  1,2*K  +  1)  =  XTEMP (C , 2 , 3 , 3 ) +XHH (C , 2  *  1 , 2*J  +  1 , 2*K+ 1 ) 
XHH(C,2*1+1,2*J+1,2*K+1)=XTEMP(C,3,3,3) +XHH (C, 2*1+1 , 2*J+1 , 2*K+1 ) 

140  CONTINUE 

DO  35  O  =  1,3 
DO  35  P  =  1,3 
DO  35  Q  =  1,3 
DO  35  S  =  1,3 

IF  (  O  .EO.  DTHEN 
DPTEMP  (P,Q,S)  =  0 
END  IF 

XTEMP (0,P,Q,S)  =  0 

35  CONTINUE 

30  CONTINUE 

DO  45  K  =  1,R 
DO  45  J  =  1,R 
DO  45  I  =  1,R 
DO  45  L  =  1,3 

IF  (XHH(L,I,J,K) .NE.O)  THEN 

XH(L,I,J,K)  =  XHH(L.I,J,K)  +  XH(L,I,J,K) 

END  IF 

45  CONTINUE 

C  COMPUTE  THE  COST 

COST  =  0 
SUM  =  0 
DO  50  K  =  1,R 
DO  50  J  =  1,R 
DO  50  I  =  1,R 
DO  50  L  =  1,3 

IF(L  .EQ.  1)  THEN 
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COST  =  COST  + 
ELSEIF  (L  .EQ.  2 
COST  =  COST  4- 
ELSE 

COST  =  COST  t 
END  IF 

SUM  =  SUM  +  XH(L 
IF(XH(L, I, J,K)  .NE 
WRITE(12, 12)L, I 
ENDIF 

50  CONTINUE 

WRITE (12. 13)  COST, SUM 
RETURN 
END 


I»XH(L, I, J.K) 
THEN 

J»XH(L, I, J.K) 
K*XH(L. I, J.K) 

I.  J.K) 

0 ) THEN 

J. K, XH (L, I , J.K) 
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